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Abstract 

Critical transitions occur in a wide variety of applications including mathematical biology, 
climate change, human physiology and economics. Therefore it is highly desirable to find early- 
warning signs. We show that it is possible to classify critical transitions by using bifurcation 
theory and normal forms in the singular limit. Based on this elementary classification, we ana- 
lyze stochastic fluctuations and calculate scaling laws of the variance of stochastic sample paths 
near critical transitions for fast subsystem bifurcations up to codimension two. The theory is 
applied to several models: the Stommel-Cessi box model for the thermohaline circulation from 
geoscience, an epidemic-spreading model on an adaptive network, an activator-inhibitor switch 
from systems biology, a predator-prey system from ecology and to the Euler buckling prob- 
lem from classical mechanics. For the Stommel-Cessi model we compare different detrending 
techniques to calculate early-warning signs. In the epidemics model we show that link densi- 
ties could be better variables for prediction than population densities. The activator-inhibitor 
switch demonstrates effects in three time-scale systems and points out that excitable cells and 
molecular units have information for subthreshold prediction. In the predator-prey model ex- 
plosive population growth near a codimension two bifurcation is investigated and we show that 
early-warnings from normal forms can be misleading in this context. In the biomechanical model 
we demonstrate that early-warning signs for buckling depend crucially on the control strategy 
near the instability which illustrates the effect of multiplicative noise. 
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1 Introduction 

A critical transition (or tipping point) is a rapid sudden change of a time-dependent system. For the 
introduction we shall rely on this intuitive notion; the mathematical development starts in Section 
[2j Typical examples of critical transitions are drastic changes in the climate \70\ [1] , in ecological 
systems [Ml El] , in medical conditions [HU [96] or in economics [50l [51] . Reviews of recent progress to 
develop early-warning signals for these critical transitions from an applied perspective can be found 
in [84:\ [83] . The goal of a mathematical theory should be to provide qualitative and quantitative 
conditions to check whether a drastic change in a dynamical system can be predicted before it 
occurs; note that it is obvious that certain transitions are very difficult to predict, for example, due 
to large noise effects [3lJ or non-smooth transitions [47] . 

A basic assumption in many applications is that the underlying process is deterministic but 
is subject to small random fluctuations. Furthermore, one often assumes that the change occurs 
rapidly in comparison to the current state dynamics. Elementary remarks how the mathematical 
theory of stochastic fast-slow systems can be used to encapsulate these hypotheses can be found in 
|65j . In particular, several one-parameter normal form models were studied and a more detailed link 
between rising variance [23], rising autocorrelation [28], time series analysis and dynamical models 
was pointed out. For additional references on critical transitions we refer the reader to ^4] and [65] . 




We outline our results without stating detailed technical assumptions. It will be assumed that 
the main dynamics near a critical transition is governed by an ordinary differential equation (ODE). 
A classification which bifurcations are critical transitions based on a definition suggested in [65] 
is explained. In a suitable singular limit this classification is a simple exercise dealing with all 
bifurcations up to codimension two. Some of the details for this classification are explained since 
one has to determine, at least once, which conditions on the fast and slow subsystems of a multi-scale 
system near higher-codimension bifurcations lead to trajectories that resemble critical transitions 
observed in applications. To model the random fluctuations stochastic differential equations (SDEs) 
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with sufficiently small white noise are used. We calculate asymptotic formulas for all possible 
covariance matrices associated to sample paths approaching a critical transition. The setup for 
the calculations is straighforward and is based on normal form assumptions, approximation by 
Ornstein-Uhlenbeck processes and the solution of a few algebraic equations. An error estimate 
for the asymptotic expansions is proven for the fold bifurcation by analyzing stochastic difference 
processes and applying elementary moment estimates, thereby avoiding more advanced techniques 
|15j for a certain regime. The focus on the fold bifurcation is justified as it is one of the most 
frequently encountered critical transitions |9H For the same reason, we also provide higher- 
order asymptotic expansions for the variance as doubly singular limit expansions with small noise 
and small time scale separation for the approach towards a fold point. 

Then we use the mathematical results in a wide variety of models. For each application the 
theoretical predictions are compared with numerical results. We briefly describe which important 
results are obtained within the examples. In a box-model of atmospheric and ocean circulation 
we test different approaches to estimate the variance from a given time series and suggest a new 
method motived by fast-slow systems. In a discrete epidemic spreading model a moment expansion 
is used to simplify an adaptive network and to analyze the onset of an epidemic. It is shown that 
predictability in adaptive networks can be improved by focusing on link dynamics instead of node 
dynamics. A model from systems biology is used to explain the effect of two critical transitions 
linked in a three-time scale systems. A predator-prey model illustrates the effect of multiple system 
parameters which can potentially hide early-warning signals that are usually expected to occur in 
ecology. The last example treats buckling of a spring in the context of a biomechanics experiment. 
The model for this experiment shows how parameter-dependent non-additive noise influences, and 
systematically changes, observed early-warning signs. The examples from epidemics, biomechanics 
and systems biology also seem to be among the flrst (or even the flrst) ones where early-warning 
signs for critical transitions have been applied in the respective flelds. 

In summary, our theoretical results combine well-known elementary mathematical tools from 
bifurcation theory, fast-slow systems, real analysis, stochastic differential equations, probability, 
asymptotic analysis, numerical continuation/integration and time series analysis to systematize 
some of the aspects of critical transitions. In this way, we make progress towards the major open 
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question to develop a unified critical transitions theory [8^. Although no complicated technical 
steps are treated we hope that our work forms a starting point to motivate new mathematical 
insights into predictability for dynamical systems; see also Section [HI Our second contribution is 
to show that abstract critical transitions theory can yield very useful conclusions with immediate 
value for applications. 



The paper is structured as follows. Section [2] describes the background from deterministic fast- 
slow systems. Section [3] contains the classification results. Section H] reviews theory for stochastic 
fast-slow systems based on which we prove the error estimate for asymptotic moment results near 
folds. In Section [5] the leading-order asymptotic scaling laws for the covariance are obtained and 
in Section E] these results are refined for the fold. Section [7] contains the five important examples. 
Section [8] provides an outlook how the framework developed here could be extended. 



Convention: Whenever a citation with detailed page numbers at the beginning of a result 
(Theorem, Lemma, etc.) is given then the statement and proof can be found in the reference. 

2 Brief Review of Fast- Slow Systems 

We recall the necessary definitions and results from multiple time scale dynamics [301 [56l [73l HQ] 
that are required to define critical transitions. A fast-slow system of (ODEs) is given by 



,f = ei; = f{x,y,e), 

ds 



— - y = 9{x,y,e), 



where < e <^ 1, x £ are fast variables and y £ M" are slow variables. The maps 
/ : M.rn+n+1 _^ ^ . _^ assumed to be sufficiently smooth. If /, g do not 

depend on e we omit the e-argument and write e.g. f{x,y) instead of f{x,y,e). Changing in ([TJ 
from the slow time s to the fast time t = s/e gives 

f = X' = f{x,y,e), 



dy ^/ 

dt 



y = (^9{x,y,e). 
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Henceforth (') will denote differentiation with respect to the slow time s and prime differentiation 
with respect to the fast time t. The singular limit e — ?• in ([T]) yields the slow subsystem 



= f{x,y,0) 



(3) 



y = 9{x,y,o), 



which is a differential-algebraic equation restricted to the critical manifold Cq := {{x,y) G 
l^m+n . f^x,y,0) = 0}. The fast subsystem is obtained as the singular limit of ^ 



where the slow variables can be viewed as parameters. The flows generated by ([3]) and ^ are called 
the slow flow and the fast flow respectively. A point p E Cq is an equilibrium point of the fast 
subsystem. The critical manifold is normally hyperbolic at p € Cq if the m x m matrix Dxf{p) 
has no eigenvalues with zero real parts. In this case, the implicit function theorem provides a map 
/iQ : — ^ that describes Cq, locally nearp, as a graph Cq = G : x = ho{y)}. Then 

the slow subsystem ^ can be written more concisely as y = g{hQ{y),y). If all eigenvalues of Dxf{p) 
are negative (positive) then Cq is attracting (repelling); other normally hyperbolic critical 
manifolds are of saddle-type. Observe that Cq is attracting at p if and only if the fast subsystem 
has a stable hyperbolic equilibrium at p. Fenichel's Theorem provides a complete description of the 
dynamics for normally hyperbolic invariant manifolds for sufficiently smooth vector fields {f,g)- To 
state the result, we recall that the Hausdorff distance between two sets V,W C M™""^" is given 



dniV, W) = max < sup inf \\v — w\\, sup inf \\v — w\\ > 
where || • || denotes the Euclidean norm. 

Theorem 2.1 (Fenichel's Theorem [35]). Suppose S = Sq is a compact normally hyperbolic 
submanifold (possibly with boundary) of the critical manifold Cq. Then for e > sufficiently small 
there exists a locally invariant manifold diffeomorphic to Sq. has a Hausdorff distance of 
0{e) from Sq and the flow on S^ converges to the slow flow as e — )• 0. S'e is normally hyperbolic 



X 



f{x,y,0) 



(4) 



by 
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and has the same stability properties with respect to the fast variables as Sq (attracting, repelling 
or saddle type). 

Sf, is called a slow manifold and is usually not unique. In regions that remain at a fixed 
distance from the boundary of S^, all manifolds satisfying Theorem 12.11 lie at a Hausdorff distance 
0{e~^^'') from each other for some K > with K = 0(1). The choice of representative will be 
irrelevant for the asymptotic analysis we are interested in here; see also ^Q\- If the choice of subset 
Sq is understood then we also write for the slow manifold associated to Co and refer to as 
"the" slow manifold. 




-0.4 -0.2 0.2 0.4 

X 

Figure 1: Critical transition at a fold bifurcation of the fast subsystem. The critical manifold Co 
splits into a repelling part (dashed grey) and an attracting part (solid black). A typical candidate 
trajectory 70 (red) is shown. We have also sketched the two different regions: (Rl, blue) where 
normal hyperbolicity of the critical manifold holds and (R2, green) near the bifurcation point. 

A candidate trajectory 70 is a concatenation of slow and fast subsystem trajectories; see 
Figured) More precisely we define a candidate as a homeomorphic image 70 (i) of a partitioned 
interval o = to < < " " " < *m = ^1 where the image of each subinterval ^o{tj-i,tj) is a trajectory 
of either the fast or the slow subsystem and the image 70(0, b) has an orientation that is consistent 
with the orientations on each subinterval 7o(tj_i,tj) induced by the fast and slow flows. Note that 
the intervals {tj,tj^i) do not necessarily correspond to the time parametrizations of the fast or 
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slow subsystem; to achieve such a parametrization one has to pick a convention such as using the 
slow time and compactifying infinite time intervals as necessary. 

If consecutive images 7o(tj-i,tj) and 7o(tj,tj+i) are trajectories for different subsystems then 
we say that 7o(ij) is a transition point. 

Definition 2.2. Let p = {xp,yp) G Co be a point where the critical manifold Co is not normally 
hyperbolic. We say that p is a critical transition if there is a candidate 70 so that 

(CI) 7o(tj--i,tj) is contained in a normally hyperbolic attracting submanifold of Cq, 

(C2) p = 70 (tj) is a transition point, 

(C3) and 7o(tj_i,tj) is oriented from 7o(tj_i) to 70 (ij)- 

Definition 12.21 was suggested in |65j. It is related to the concept of "hard" or "catastrophic" 
loss of stability ([69], p. 87 or [5], p. 36) but does not coincide with it. Note that Definition 12.21 is 
entirely based upon the singular limit e = 0. Definition 12.21 is simple, easy to verify for a system, 
concretely includes the focus on the candidate orbit occuring in an actual time series and also seems 
to represent all the requirements laid out in [84J. Note carefully that (CI) excludes slow canard 
orbit segments in repelling parts of the critical manifold but see Section [8] for possible extensions. 

Proposition 2.3. Suppose p = (xp, Up) is Lyapunov stable with respect to the fast subsystem, then 
there is no critical transition at p. 

Proof. Suppose 70 is a candidate that satisfies (CI) of Definition 12.21 If p = jo{tj) is a transition 
point then the orientation condition (C3) implies that the 7o(tj,tj+i) is an orbit segment in the 
fast subsystem starting from p. Since p is Lyapunov stable we have reached a contradiction. □ 

In this paper, we are interested in the approach of trajectories to critical transitions as illustrated 
in Figured! This approach towards a critical transitions can be subdivided into two main regions: 
(Rl) Fenichel's Theorem applies near a normally hyperbolic critical manifold and (R2) Fenichel's 
Theorem fails near the bifurcation point. Note that Fenichel's Theorem implies that near a local 
bifurcation point (x, y) = {xp, yp) of the fast subsystem the region (R2) shrinks to (xp, yp) as e — )• 0. 
By making e sufficiently small, we should start to focus on (Rl). Whenever we consider decreasing 
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y — )• yp we shall make the assumption that e has been chosen small enough so that we stay inside 
region (Rl). For example, for a fold point ^64j it is known that the size of (R2) scales like 

(x,y)~(0(ei/3),0(e2/3))eM2 

as e — )• 0; see also Lemma 16.11 Therefore we would assume that e^/^ <^ y as y gets small. We 
formalize our assumptions by restricting the analysis to a compact domain contained in the region 
(Rl): 

(AO) Fast-slow systems will be considered on a compact domain for 2?(e) = T> = x Dy c M"* x 

that depends smoothly on e. T>{e) is chosen so that an attracting slow manifold is contained 
in I'(e), the intersection dD{e) n is transverse and T){e) is contained in the basin of 
attraction of Cf; the slow manifold Cf is given locally as a graph = {{x^y) G ^(e) : 
X = he{y)}, for /i^ : Dy — t- D^- Furthermore, fast subsystem local bifurcation points will lie 
on 9P(0) and asymptotics with respect to y — )• yp is chosen depending on e so that normal 
hyperbolicity holds; see Figured! 

Note that this means that all scaling estimates we derive are restricted to a bounded domain. 
Within this bounded domain no other attracting critical manifold perturbs to a slow manifold. 

3 Fast Subsystem Normal Forms and Critical Transitions 

We assume familiarity with the normal form approach to bifurcation theory ( |43j . p. 138) and apply 
it in the singular limit to the fast subsystem viewing y E M" as parameters. The number of 
slow variables y G M*^ is chosen as the codimension of the bifurcation. We are going to check 
which bifurcations are critical transitions in the sense of Definition 12.21 Note carefully that this 
classification, although complete on the singular limit level e = 0, has interesting possible extensions 
which are discussed in Section [8l 

Assume without loss of generality that the bifurcation point is at x = (0, . . . , 0) =: and 
y = (0, . . . , 0) =: 0. To reduce the analysis to normal forms we will assume that all the necessary 
genericity conditions (non-degeneracy and transversality) are satisfied |69j : 



8 



(Al) / G C^(M™'+"+^, M*") where r > 1 is chosen according to the differentiabihty required by 
normal form theory. We also assume that g G C^(M™'^""*"^, M"). 

(A2) The genericity conditions for bifurcations hold so that normal form theory applies. 

The only generic codimension one bifurcation with x G is the fold bifurcation with 
normal form ([69j, p. 84) 

f{x,y) = -v-x^. (5) 
Considering the dynamics of the slow variable as g{x, y) [65] we get the fast-slow system 

/ 2 

X = —y — X , 

(6) 

y' = (^9{x,y), 

where y G since we have a codimension one bifurcation. The critical manifold for ([6]) is Co = 
{{x,y) £ M? : X = zt-^/— y =: h^{y),y < 0}. Cq Pi {x = \/—y, y < 0} is attracting and Cq n {x = 
—\J—y^ y < 0} is repelling. The linearization Dxf{hQ{y),y) around the attracting branch of the 
critical manifold will be crucial for the calculations in Section [5] as it describes the dynamics in the 
region (Rl); therefore we will calculate/record this linearization for each critical transition. 

Lemma 3.1. If g{0,0) > then ^ has a critical transition at {x,y) = (0,0). 

The proof is obvious and similar results hold for the pitchfork and transcritical normal forms 

f{x,y) = yx + sx^ Dxf{ho{y),y) = y, for s = ±1, p.284), (7) 

fix,y) = yx-x\ D,f{ho{y),y) = y. ([43], p. 149)). (8) 

Lemma 3.2. If g{0, 0) > then ([7]) has a critical transition at (x, y) = (0, 0) if and only if the 
pitchfork bifurcations is suhcritical (s = 1). If g{0,0) ^ then ([8]) has a critical transition at 
ix,y) = (0,0). 

The remaining one-dimensional fast subsystem is the codimension two cusp bifurcation. The 
normal form is ([69], p. 304-305) 

f{x,y) = yi+y2X + sx^, for s = ±1 (9) 
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where the fast dynamics x' = f{x, y) is augmented with two-dimensional slow dynamics y' = 
eg{x, y), y = (^1,^2) G I^^- The critical manifold for ([9]) is Co = {(x, y) G M'^ : = yi + y2X + sx'^}. 
Due to the two-dimensional slow flow it is slightly less obvious to determine under which conditions 
the cusp bifurcation is a critical transition. 

Lemma 3.3. There is no critical transition for Q at {x,y) = (0,0) if s = —1. If s = 1 then ([9]) 
has a critical transition at {x,y) = (0,0) if and only if §2(0,0) > and gi{0,0) = 0. 

Proof. First consider the case s = — 1. At yi = = ^2 the fast subsystem is x' = —x^. It is easy to 
see that x = is asymptotically stable and Proposition 12.31 implies that there cannot be a critical 
transition at (x,y) = (0,0). For s = 1 the fast subsystem is x' = x^ so that a candidate orbit 70 
can have a segment 7o(tj,tj+i) in the fast subsystem oriented from 70 (tj) to 7o(tj+i). To see that 
there exists an attracting critical manifold connecting to (x, y) = (0, 0) we need the unfolding of 
a cusp bifurcation. The linearization of ([9]) is Dxf{x,y) = y2 + 3x^ and the stability of the slow 
manifold changes at fold points when Dxf\co ~ ^- Given the two equations 

= 2/1+ y2X + x^ 

= 1/2+3x2 

the variable x can be eliminated which yields the classical cusp curve. After a projection into 
the (yi, y2)-plane it is given by F := {(2/1,2/2) G '■ ^yf + '^"^Ui = 0}- ^ repelling subset of the 
critical manifold is := Cq n {4y| + 27yf > 0}. The set Cq n {4yf + 27yj < 0} splits into three 
branches corresponding to the three solutions of f{x,y) = where two branches Cq^ are repelling 
and one branch Cq is attracting; observe that 2/2 < for any of the three branches. Now consider 
a candidate 70 with 7o(tj-i,tj) C Cq = {x = /io(2/)}- The slow flow on Cq is given by 

2/1 = 51(^0(2/), 2/), ^^^^ 
2/2 = g2{ho{y),y). 

Since the 2/2-component of 7o(tj-i) is negative we must have 2/2(0,0) > 0. Furthermore, we know 
that trajectories of (fTOj) on reach (2/1,2/2) = (0,0) if and only if 4yf + 27yl < holds for the 
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y-components of 7o(tj_i, tj). Considering the branches of T we have 

In particular, we find that hmjy2^or±(y2) = which imphes that gi(0,0) = if the candidate 70 
reaches {x, y) = (0, 0). Hence there exists 70 as required by Definition 12.21 if and only if (72(0, 0) > 
and 51(0,0) = 0. □ 

The attracting part of the critical manifold Cg (as introduced in the previous proof) is Cg = 
{(x, 2/1,1/2) G : a; = /ig(y)}. The linearization is 

D,f{ho{y),y) = y2 + 3ho{yf = Oy{y2) as y ^ 0. (11) 

where the notation Oy{-) indicates asymptotic scaling as y — )• under the assumption that Fenichel 
Theory is still valid; see also Section [2] where this is referred to as region (Rl). The asymptotic 
scaling in ([TT]) holds since points on Cg satisfy 1/2 + < because on Cg we have 

y2 < 0, x€ [-(-y2/3)'/^ (-^2/3)^/'], yi = -xy2 - x^ 

Therefore, x = ho{y) grows at most like ^/y2 as y — t- and the scaling law in (jlip follows. This 
concludes our discussion of one-dimensional fast subsystem bifurcations. 



For two fast subsystem variables consider the codimension-one Hopf bifurcation normal form 
([69], P.98) 

fi{x,y) = yxi-X2 + hxi{xl + xl), 

{12) 

/2(x,y) = xi+yx2 + hx2{xl + xl), 

where ^i is the first Lyapunov coefficient. The critical manifold for (|12p is Cg = {(x,y) S 
: x = 0} where Cq H {y < 0} is attracting and Co H {y > 0} is repelling and the linearization is 



D,f{0,y) 



y -1 
1 y 



(13) 
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Lemma 3.4. // g{0, 0) > then (jl2p has a critical transition at (x, y) = (0, 0) if and only if the 
Hopf bifurcation is subcritical (li > 0). 

For vanishing first first Lyapunov coefficient {li = 0) a codimension-two generalized Hopf (or 
Bautin) bifurcation occurs with normal form ([69], p. 313) 

fi{x,y) = yixi - X2 + y2Xi{xl + xl) + l2Xi{xl + xlf , 

(14) 

f2{x,y) = xi+yiX2 + y2X2{xi + xl) + l2X2{xl + xiy, 

where I2 = ±1 is the second Lyapunov coefficient. The critical manifold is Cq = {{x,y) G : 
X = =: holy)}. The linearization Dxf{hQ(y),y) coincides with the linearization (|13p for the Hopf 
bifurcation upon replacing y by yi. 

Lemma 3.5. The Bautin bifurcation is not a critical transition if I2 < 0. 

Proof. Without loss of generality let I2 = —1 then the fast subsystem at (1/1,2/2) = (0)0) is 

x[ = -xi{xl + xlf =: fi{x), 

(15) 

x'2 = -X2{xl + xlf=:f2{x). 

where / = (/i, /2). Define a function y : — )■ M by V{xi,X2) := + x^- Observe that V{x) > 
for X 7^ and 

jV{x) = D^V{f{x)) = -{2x1 + '^xDixl + xlf < 

for X 7^ 0. Therefore V{x) is a Lyapunov function and x = is asymptotically stable as an 
equilibrium point of (|15p and Proposition 12.31 finishes the proof. □ 

Lemma 3.6. If I2 > then the Bautin bifurcation is a critical transition if and only if gi{0, 0) > 
and either (a) 52(0,0) or (b) 52(0,0) = and §§(0,0) < i. 



Proof. The critical manifold splits into two 2-dimensional planes Cg = Co H {yi < 0} and Cq = 
Co n {yi > 0} where Cg is attracting and Cq is repelling. The condition 51(0,0) > implies that 
the slow flow y = 5(0, y) has trajectories that start in Cq and reach (x, y) = (0, 0) in finite time. 
This guarantees the existence of a candidate 70 satisfying (CI) and (C3) of Definition 12.21 The 
proof of Lemma 13.51 shows . upon reversal of time in equation (|15p . that x = is an asymptotically 
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unstable equilibrium point of the fast subsystem at y = when I2 = 1- We note from the unfolding 
of a Bautin bifurcation (see [69], p. 314) that saddle-node bifurcations of limit cycles for the fast 
subsystem occur on the curve LPC := {(2/1,^2) G : juf = Hi, 1/2 < 0}. The conditions on the 
slow flow guarantee that the candidate orbit enters the fast subsystem region without limit cycles 
and with an unstable equilibrium point; this region is given by 

|(yi,y2) G M2 : y2 < 0,yi > ^y?| U {(yi,y2) G : y2 > 0, yi > O} . □ 

The last codimension two bifurcation with two fast variables is the Bogdanov-Takens bifurcation 
with normal form fj43j. p. 365) 

fi{x,y) = X2, 

(16) 

f2{x,y) = yi + y2X2 + xl + SX1X2, 

where s = ±1. The critical manifold is Co = {(x, y) G : X2 = 0, xi = it-^— yi} so that we always 
require yi < 0. 



Lemma 3.7. The Bogdanov-Takens bifurcation (jl6p is a critical transition for s = —1 if and only 
if 92(0,0) > and gi{0,0) = and for s = 1 if and only if (a) gi(0, 0) > or (h) gi(0, 0) = 0, 
52(0,0) >0, ||(0,0) < -2 . 

Proof. As usual we consider the fast subsystem at y = 

x\ = X2, 

(17) 

x'2 = xf + SXiX2- 

The theory of non-hyperbolic equilibria in planar analytic vector fields (|80j, p. 151; see also [2]) 
implies that x = is a cusp point. Hence there exists a candidate 70 with a fast-subsystem orbit 
segment 7o(tj,tj+i) oriented from 7o(ij) = (0,0) to 7o(tj+i). It remains to show when we can 
approach (x, y) = (0, 0) via the slow flow on an attracting critical manifold. The linearization 
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around the critical manifold is 







1 



(18) 



V 




with Tr(L)x-/|co) = Sy/—yi and det(I?2,./|co) = — 2(ib-^— yi). Hence the critical manifold is 
attracting if and only if ?/2 i SyJ—yi < and 2{±^/—yl) < 0. For s = —1 this yields the set 

Co = {(.xi,X2,yi,y2) G : yi < 0, 2:2 = 0,xi = -^f^,y2 < -V^}- 

A candidate orbit starting in Cq will reach (x, y) = (0, 0) if and only if 52(0, 0) > and 51 (0, 0) = 
where the second condition is required since the curve {2/2 = ~\/~yi} approaches the origin 
tangentially i.e. ^ (—7/2)1^=0 = 0. The second case for s = 1 is similar and follows using the 
symmetry {xi,X2,yi,y2,t) ^ {xi, -X2,yi, -y2, -t). □ 

The linearization for the Bogdanov-Takens bifurcation has already been recorded in (jlSp but 
notice that we must have the condition y2 < ^s^/—yl to be on Cq. The asymptotic result as y — )• 
for the linearization is 



The two remaining codimension-two bifurcations (fold-Hopf and Hopf-Hopf) require three and four 
fast dimensions. The fold-Hopf (or Gavrilov-Guckenheimer) bifurcation has normal form (|69j. 



fi{x,y) = yi + xl + s{xl + xl), 

f2{x,y) = y2X2 - UJX3 + 9xiX2 - X1X3 + xlx2, 

fz{x,y) = ujx2 + y2X3 + xiX2 + 0x1x3 + x\x3, 
where s = ±1 and 6 = 9{y) satisfies 0(0) ^ and a; 7^ 0. The critical manifold is given by 




(19) 



p.338) 



Cq = {(x, y) ■.X2 = Q = X3, XI = ±^/^, yi < 0}. 



Lemma 3.8. The fold-Hopf bifurcation is a critical transition 
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• for e{0) > 0, s = 1 if and only if (a) gi{0,0) > or (h) 51 (0,0) = and 52(0,0) > 0, 

• for 6'(0) > 0, s = -1 i/ and only if gi(0, 0) > and 52(0,0) < J2(0,0) where ^2(0,0) is the 
y2-component of the tangent vector to the "cycle blow-up curve" (cf. 169^ . p. 343), 

• for 61(0) < 0, s = 1 if and only i/ 51 (0,0) = and 52(0,0) > 0, 

• for 61(0) < 0, s = -1 if and only i/ 51 (0,0) = and 52(0,0) > 0. 

Proof. The same techniques as before will apply so we just sketch the proof. The fast subsystem 
at y = is 



xl + s{xl + xl), 



^'3 



2 = -UJX3 + 9xiX2 - X1X3 + xlx2, 

0JX2 + X1X2 + 6x1x2, + xfxs, 



(21) 



Changing to cylindrical coordinates (x2,X3) = (r cos r sin </>) in (|'2ip and neglecting the angular 
component (j), since it is always a neutral direction with respect to attraction and repulsion for the 
critical manifold of equilibrium points, we get a two-dimensional system 



xf + sr , 
r{6xi + x1] 



(22) 



It can be checked that the origin (xi,r) = (0,0) is unstable for (j2ip . Therefore we can find a 
candidate that leaves the bifurcation point in a fast direction. The attracting part of the critical 
manifold is computed from the linearization 



D.flco 



I 



\ 



±2^^ 

52 ± 6V-yi T V-yi 

uj± 52 ± J 



(23) 



and is given by Cg = Co n {x = —\/—yi, 52 < ^\/— 5i}- The conditions on the slow flow y = g = 
(91,92) can be derived from the unfolding of the fold-Hopf bifurcation (see [69], p. 339-345). □ 



The linearization is given by (j23p : we note that the condition of the approach via Cn means 
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that the leadmg order approximation to Dxf\concg as (yi,?/2) — ^ (0 ,0) is given by 



Dxflconcg 



±2^^ 

Oyi^/^i) -uj 
u Oy{^^) I 



(24) 



As the last case we are going to consider is the Hopf-Hopf bifurcation. We shah not discuss 
the complicated unfolding ([69j, p. 351-370; |43) . p. 396-411) in detail to show when the Hopf-Hopf 
bifurcation is a critical transition. A normal form in polar coordinates (j"!, ^2, ^i, ^2) = ^) is 
([69], P.358) 

fi{r,0,y) = ri{yi +purj +pi2r^ + sir^), 



(25) 



Mr,0,y) = wi, 
f4{r,9,y) = UJ2, 

where 101^2 are the imaginary parts of the eigenvalues at the bifurcation point y = 0; pij and si^2 are 
further parameters. We note that all parameters also depend on the slow variables y but are usually 
assumed to be non-zero at the bifurcation point. Observe from (I25p that the critical manifold is 



Co = {{r, 0, y) e X m2 . = = rs} = {{x, y) e 



: Xj = for j = 1,2,3,4}. 



To study whether candidate orbits can leave the fast subsystem for y = we would have to study 
the nonlinear stability of the origin depending on the parameters. It is not difficult to see that the 
Hopf-Hopf bifurcation is not always a critical transition depending on parameter values but there 
are cases when it is a critical transition. Instead of providing this detailed study (which can be 
inferred from the unfoldings in |69] ) we shall only state one important linearization 



Dxf{x,y)\conC^ 



^ yi -uji ^ ^ 

wi yi 

y2 -^2 

\ UJ2 y2 J 



(26) 
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It will always be assumed for (j26p that u)i^2 7^ 0; note that we also exclude resonances koji = luj2 for 
A; + / < 3 as non-resonance conditions are non-degeneracy conditions for the Hopf-Hopf bifurcation 
(cf. assumption (A2)). We remark that the main bifurcation phenomena of interest near a Hopf- 
Hopf bifurcation are global orbits (limit cycles and tori) which would not be captured by our local 
analysis anyway. 



Name 


N.-Form 


Critical, if... 


Dxflconcii =■ Ao{y) 


Fold 


Eq. dSD 






Pitchfork 


Eq. dZD 


s = 1, 5 > 


y 




Fn 




y 


Hopf 


Eq. (m 


li>0,g>0 




Cusp 


Eq. dSD 


s = -1, 51 = 0, 52 > 


Oyiy2) 


Bautin 


Eq. dm 


''2 ^ u, yi ^ u 

and (a) 52 / or 

(b) 52 = 0, dy,g2 < 1/2 




Bog.-Tak. 


Eq. (HI]) 


s = -1, 51 > 0, 52 = 
s = 1 and (a) 51 > or 
(b) 51 = 0,52 >0,a^, < -2 


{ 


1 ^^ 


Fold-Hopf 


Eq. ^ 


< 0, 51 > 0, 52 = 

9 > 0, s = 1 and (a) 51 > 0, 

or (b) 51 = 0, 52 > 

(9 > 0,s = -1,51 > 0,52 < J2 


/ ±2V-yi \ 
0,(V-5i) 


Hopf-Hopf 


Eq. ^ 


special case only 




/ 51 -wi \ 
yi 
52 — a;2 

\ W2 52 / 





Table 1: Results for fast subsystem bifurcations. The additional hypotheses on the slow flow 
y = g{x,y) at {x,y) = (0,0) are abbreviated and we always understand gj as 5j(0,0) for j = 1,2 
and 5 as 5(0, 0) in this table. The Hopf-Hopf bifurcation has not been analyzed in detail and only 
a particular case is stated. The last column records the linearization around the attracting branch 
of the critical manifold. 



Having finished the exercise it is now clear which local fast subsystem bifurcation points are 
critical transitions under suitable slow flow conditions. We record the results developed in Lem- 
mas [3?1][3]8] as well as the resulting linearizations -Dx/lconC" in Table [1] where we introduced the 
shorthand notation ^0(5) := -Da/(/io(y), 5) = D^flconcg- 

Let us point out again that the classification results are for the singular limit e = 0. Detailed 
unfoldings for the deterministic case for e > are known for the fold, pitchfork, transcritical and 
Hopf bifurcations [64:\ 162] 163] [76] . Partial results are available for the Bogdanov-Takens bifurcation 
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|25j and the cusp [22] is work in progress. Section [8] provides an overview where future work is 
needed. 



4 Sample Paths and Moments for Stochastic Fast-Slow Systems 

Let {Ws}s>o be a /c-dimensional Brownian motion on a probabihty space {0,,J-,¥). Consider the 
fast-slow stochastic differential equation (fast-slow SDE) 



dxs = jf{xs,ys)ds + ^F{xs,ys)dWs 
dys = g{xs,ys)ds. 



(27) 



which is understood as an Ito-SDE |79j . Noise acting on the slow variables y will not be considered 
explicitly but it is implicitly included in all of our results as it appears as a higher-order term ( [TS] , 
p. 145; [65], p. 1026). In addition to the assumptions (AO)-(Al) that hold for the deterministic part 
of (j27p we require the following hypothesis: 



(A3) F £ C72(j^m+, 



n Tn>mxk\ 



and the noise level a = o"(e) depends continuously on e. 



(A4) We consider small noise with \im^^oa{€) = 0. 

To understand the effect of a deterministic smooth invertible normal form transformation (coor- 
dinate change) u{x, y) = {X, y) G M™ x M", with u G ^"(R™ x M", R"^ x M") we need the following 
result which directly follows from Ito's formula ([79j, p. 44). 

Lemma 4.1. Consider the fast variable equation for ([2] 



dxs = -f{xs, ys)ds + -^F(xs, ys)dWs 



then, using the notations Zg = {xs,ys) o-nd Zg = {Xg^Ys), we have 



dXi'^ 



^/«(X„y,) + 0(1) + O ds + ■^^F^'HXs,Ys)dWs 



ds + ^F^'\u-\Zs))dWs 



(28) 



where superscripts (i),(i) denotes the i-th resp. j-th row/component. 
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Since g{xs,ys) = is smaller than the first term, the only term that could be of leading 

order and obstruct the transformation to normal form for the deterministic part of (I27p is of order 
0{a'^/e). By (A4) we have cT^(e)/e ^ 1/e which implies that the third term is also of higher- 
order after the normal form transformation in comparison to the deterministic 0(l/e)-term. We 
now formally truncate ()28p by discarding the two terms of order lower than 0{l/e) as well as the 
polynomial terms appearing in f^^\Xs,Ys) which are of higher-order than the leading normal form 
terms (e.g. ior the Md -Y + 0{Y^, XY, X^) + 0{e) + 0{a^ /e) ^ -Y). On the basis of this 
formal truncation we now work with (1270 where / is chosen from the set of deterministic normal 
forms discussed in Section [3l There are several interesting remarks regarding the formal truncation; 
see also Section 8. 

Remark 1: In the deterministic case on the fast time scale, discarding higer-order polynomial 
and 0{e) terms is well understood for generic fast subsystem bifurcations as shown e.g. in ([90], 
Proposition 2.1, Section 4.1; |62j, equation (2.5), section 2.4). Intuitively this is also clear since 
small perturbations do not change the unfolding of 1- or 2-parameter generic bifurcations for gen- 
eral smooth vector fields [99j . 

Remark 2: For the deterministic (cr = 0) pitchfork and transcritical bifurcations, which are 
not generic for general smooth vector fields, the C'(l)-term in ()28p is relevant as shown e.g. in 
(|63j. Lemma 2.1, Theorem 2.1) in a region of the type (R2) near the singularity. The stochastic 
early-warning signs in a normally hyperbolic attracting region, such as (Rl), are not expected to 
depend upon these terms (see the discussion of the attracting regime in jT2]) but we do not verify 
this here and work with the formal truncation. 

Remark 3: It might be possible to weaken the assumption (A4) and to give a rigorous proof for 
a suitable 'equivalence' of a given SDE and its normal form. For the deterministic case, topolog- 
ical equivalence is known but for the stochastic case one needs different concepts such as random 
normal form transformations as suggested by Arnold and co-workers [3]. 
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Once the SDE (j27p has been transformed into normal form we study sample paths {xs,ys) 
that solve (f27|) as suggested in [H]. By (AO) there exists a deterministic attracting slow manifold 

= {(x,y) G X : X = h,{y)} 

for : Dy — )• D^- The deviation of sample paths from this deterministic slow manifold is = 
Xs — hf{ys) and the variational SDE for is 



d^s = dxs - Dyh^{ys) dys (29) 
= - [f{he{ys) + (s,ys) - eDyh^{ys) g{K{ys) + (s,ys)] ds + ■^F{h,{ys) +is,ys)dWs. 



Let y^*^* denote the deterministic solution of p7p (i.e. a solution for o" = 0). For .^^ = the drift 
term in ()29p satisfies the invariance equation [lOOj for a slow manifold 

f{K{yf'), yf) - eDyK{yf') g{K{yf'),yf') = 0. (30) 

Linearizing (i29]l around = and replacing yg by y^*^* yields a lowest-order system for the process 
(.S.i,ys) given by 

d^i = \[D.f{K{yf'),yf')-eDyK{yf') D^g{K{yf'),yf')]i[ds 

+^F{K(yf').yf')dW,, (31) 
dyf' = g{h,{yf'),yf')ds. 

For notational simplicity we let 

A,{yf') := D,f{K{yf'),yf')-eDyK{yT)D^a{K{yf'),yf'), (32) 
F,{yf') := F{K{yT),yf'). (33) 



Note carefully that for e = we get the matrix Ao(y) = Dxf{hQ{y),y) which is precisely the 
linearization recorded in Table [TJ We shall always assume that initial conditions for (j3ip are 
deterministic and given by (^q, yo) = (0, yo) which corresponds to starting on the deterministic slow 
manifold. Now we can state an important result about the covariance Cov(^^) of the linearized 
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process. 

Lemma 4.2 ([15], p. 146-147). Let X., := a''^Cov{i[) then Xg satisfies a fast-slow ODE 

eX = My)X + XA{yf+F,{y)F,{yf, 

(34) 

y = 9{he{y),y)- 

Furthermore, for < e ^ 1, the critical manifold for (j34p is attracting for y S T>y and given by 

Co = {X€ M.^''^ : Ao{y)X + XAo{yf + Fo{y)Fo{yf = 0}. 

Fenichel's Theorem provides an associated slow manifold Ce = {X = He{y)} for : M" — t- 
j^mxm_ Assuming that the matrix H^(y) is invertible and that the operator norm ||ff^^^(y)|| is 
uniformly bounded for y S Dy one can define the covariance neighborhood 

B{r) := {{x,y)eV:[x- K{y)f ■ H,{y)-\x - K{y)] < r^} . 

Define the first-exit time of the original process {xs,ys), starting at s = sq, from a set A as 

TA := inf{s G [sq, oo) : {xs,ys) i A, (xq, yo) G A^} 

where A is chosen so that is a stopping time wrt the filtration generated by {(x^, ys)}s>so- 

Theorem 4.3 ([l3], p. 149-150). Sample paths stay inside B{r) with high probability. More precisely, 
there exists K{s,e,a) and k > such that IP{Tg(-j.) < min(s,rx)y)} < K{s,e,a)e~'^^^^^'^'^^\ where 
the pre-f actor K{s, e, a) grows at most polynomially in its arguments as (e, a) — )■ (0, 0) and s — t- oo. 

The main conclusion of Theorem 14.31 is that sample paths near normally hyperbolic attracting 
slow manifolds are met ast able i.e. they stay near the manifold for exponentially long times except 
when the slow dynamics moves the system near a fast subsystem bifurcation point so that the 
stopping time Txiy is reached. Theorem l4.3l does not immediately guarantee that we can use moments 
from the linearized process (}g to approximate the moments of the nonlinear process ^s. For a 
approach to this problem re-consider the general fast-slow SDE (j27p . The associated slow flow 
ODE is dy^g = g{ho{yg),y^)ds. Define x^ := ho{y^) and observe that the solutions {xs,ys) of ([27|) 
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depend implicitly on e. A complementary result to Theorem l4.3l bv Kabanov and Pergamenshchikov 
provides a convenient convergence in probability of the process {xs,ys) to (x^,?/^) as e — )• 0. 

Theorem 4.4 ([ST], p. 45-46). Suppose (A0)-(A2) and (A4) hold. We start at s = sq and consider 
a final time S > such that (xg, Us) has not left V. Then for any s E [sq, S] 

I Ol P r, J I Ol P n 

sup \Xs — Xg\^\j and sup ly^ — J-U 

0<s<5 0<s<S' 

p 

as e — 7- where — )• indicates convergence in probability. 

As a direct corollary to this result the linearized process also approximates in probability 
as both processes tend to the same deterministic limit as e — )• 0. 

/ P 

Proposition 4.5. Under the assumptions (A0)-(A2) and (A4) we have supo<s<5lCs ~ Cs\ ~^ 0; 
e — 7- 0. In particular, we have convergence in distribution as e — )• 0. 

Proof. Observe that Theorem (j4.4p can also be applied to the processes and instead of Xs 
with /lo(ys) = = 0- This yields 

sup ies-eii= sup + < sup \Cs-o\+ sup iei-o|^o. □ 

0<s<S 0<s<S 0<s<S 0<s<S 

Proposition 14.51 only states that the two stochastic processes converge to the same deterministic 
process as e — t- 0. However, for a metastable approximation one must check how the A;-th moment 
approximation depends on the time s and the time scale separation e. In particular, we are 
interested in the first and second moments and let 

5i(s,e) :=E[6]-E[ei], <52(s,e) := Cov(6) - Cov(ei). 

It is certainly possible to adapt previous results such as the work by Berglund and Gentz [15] to 
achieve explicit moment bounds. However, the techniques are rather complicated and based upon 
martingale methods, Bernstein-type inequalities, subdivison of suitable time intervals and calcu- 
lating new explicit moment bounds and aim to control probabilities path-wise. Here we are going 
to develop a very short and essentially 'algorithmic' argument for moment bounds for truncated 
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normal forms. The technique is elementary and only uses a suitable difference process, well-known 
even-moment bounds and the Cauchy-Schwartz inequality; this approach may even have the po- 
tential to simplify calculations for path-wise control such as ([lOj, Section 4). 

We shall only discuss moment approximation for the fold bifurcation which provides an outline 
how moments can be controlled in the general case. The simplest normal form model for a fold 
bifurcation with additive noise is 

dxs = -(-Vs - x'i)ds + -^dWs, 
dys = 1 ds, 

where we can also view = (s — sq) + Vso as a time variable. The attracting critical manifold 
is Cq = G : X = ^J—y = ho{y)} with an associated slow manifold C° = {x = he{y) = 

^o(y) + C'(e)}. Note that for ([35]) we have ys = Vg^^- Therefore we get that (p9]) is given by 

d^s = ^i-2^,^s - e + Oie))ds + ^^dWs, ^ ^ 

^ (36) 

dys = 1 ds. 

where we are going to formally drop the higher-order term C'(e)-term from now on. The linearized 
problem ([3T]) is 

d^s = ^i-2V^s)esds + fJWs, ^^^^ 
dys = 1 ds. 



To analyze the transient behavior we consider the difference process Vs '■= £,s — Ci- It satisfies 
the differential equation 

dvs = ^[-Vs - KiVsf - 2he{ys)vs - Cs - eDyh^{ys)]ds 

(38) 

= \[-2^sVs-es+0{e)]ds. 

where the 0(e)-term will again be dropped. We always consider a initial condition yo at time so = 
such that ys = s + yQ and < for s € [0, s*] for some s* > such that ys remains in the compact 
region T>. 



23 



Lemma 4.6. The expected value of the difference process satisfies the ODE 

^E[vs] = ^ (-2V=yI nv,] - Eif,]) . (39) 

Proof. Substract (f37|l from ([361) and take the expectation. □ 
By the variation of constants formula ([45], p. 82) the solution of (139p is given by 

E[vs] = E[vo] X{s, so) - r ^X{s, r)dr (40) 

J so £ 

where X{s, r) is the fundamental solution of ^E[ws] = ^{—2^J—ys E[vs\). If we can show that E[^1\ 
is "small" then (j40p provides a way to show that the mean of Vs remains small as well. 

Lemma 4.7 ([57J, p. 20-25). Suppose the stochastic differential equation dXg = a{Xs, s)ds + 
l3{s)dWs with X G M"* and /3{s) G M™-x^ satisfies for s G [si,S2] the stability condition X'^ a{X , s) < 
— and has uniformly bounded noise term sup^gj^^ ,,2] ||/5(s)|| < M then 

E[X'/]<pl^—j . (41) 

Applying Lemma HTTl to = Xg and equation ([36l) we see that k = 0{l/e) and M = \/a/e. 
Therefore using (I4ip with p = 1 yields 

m,p(,,^o[^)=o{^) (42) 

where s G [0, s*] to assure normal hyperbolicity with k = ©(l/e). Using the estimate (I42p in ()40p 
and assuming that E[vq\ = E[^o — ■^ol = we get the inequality 



\E\v,]\ < 







dr = (5i(s,e). (43) 



In particular, we can use the linearized process to approximate the mean 

|E[t;,]| = ms]-E[^i]\<6i{s,e). 
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Next, we define Vg := Var(,^s) — Var(,^^). 



Lemma 4.8. The difference process Vs for the variance satisfies the ODE 



(44) 



Proof. A direct calculation using Ito's formula ([86], p. 87) shows that 



^Var(6) = 2E[l(-2V^es-a(6-E[6])] 
^Var(a = 2E[i(-2V^ei)(Ci-E[^i])]+^. 



(45) 



Then using Y&t{^s) = H^'i] - H^s? and Var(^^) = E[{^if] - E[^^]2 gives 



□ 



Lemma 4.9. \E[CtCs]\ < 0{a^) and \E[Ci]E[Cs]\ < 0{a^) 



Proof. By a combination of Lemma (14. 7p and the Cauchy-Schwarz inequality its follows that 



The second results is proven similarly. 



□ 



As in the derivation of the bound (1431) we now use Lemma \4M and Lemma [4^ to conclude that 



\Var{Cs)-Var{es)\ = \Vs\ < 



0{^]X(s,r] 



dr = 62(3, e) 



(46) 



where X{s,r) is the fundamental solution of ^E[ys] = ^{—4y/—ys E[Vs]). The estimates (|33]) and 
16]) require the fundamental solutions of systems of the form 



—w(s,r) = y/—s — Uo w(s,r), w(r,r) = 1 

ds e 



w{s, r) = [i~^-yo?''-i-r-y,f/^] _ (47) 



where k > is a constant; here k = 2, 4 for the first and second moment estimates. We remark 
that in (j47p the formal condition (— s — yo)^^"^ ~ e with s = yields the critical scaling y ~ e"^^^ 
as expected from the loss of normal hyperbolicity near a fold ([62], p.291; [15], p. 87). To estimate 
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6i{s, e) and 52{s, e) one must consider the integral 




:^(^)/^(ir with ip{r) : 



2k r 



3 L 



{-s - yof' - i-r - y,f'' 



(48) 



Note carefully that (j48p has asymptotics that can be determined via Laplace's method (see [9j, 
p. 265-267). If no formal truncation, e.g. in ([36|) . is applied there are much more detailed results 
available in jlO] using explicit calculations where Laplace-type integrals still appear [15]. However, 
it seems that the ideas used here utilizing the difference process, the direct moment estimates from 
Lemma 14.71 and the Cauchy-Schwarz inequality are a simple, and essentially algorithmic, shortcut 
to lead to a Laplace-type integral. 

Proposition 4.10. Suppose {-s + yo) = ©(e^") with q < 1/3 then e'^^^'^/'dr ~ e^~" os e ^ 0. 

Proof. One calculates that (p'{r) > for r E [0, s] if s < s*. Then applying the standard asymptotic 
Laplace approximation at the endpoint s (see |9j, p. 266, (6.4.19b)) yields 



Hence we obtain from (j43|) . (I46p and Proposition 14.101 that in the normally hyperbolic regime 
V with y = 0(e^") and a < 1/3 the moment estimates are 



Lemma 14.21 gives for the fold bifurcation the desired moment approximation for the linearized 
process Var(^f,) = (T'^[H^{y)\ so that the approximation result for the variance is 



For a = the process is at an 0(l)-distance from the critical transition point at the fold and 
Var(^s) = (T(e)^[ffe(y)] + O (cr(e)^). As expected, the estimate of variance becomes less accurate 
the closer sample paths move towards {xp,yp) = (0,0). For a > 0, the error term in formula (j49p 
is asymptotic if and only if ^ a^e'" or e" ^ a. 




as e — 0. 



□ 



6i{s,e) = Oia^e-") and 52{s,e) = ©((T^e""). 




as e — 0. 



(49) 
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Note that if e'=o+"i/fc(,(y) = 0{a) for all A; > /cq > then 



Var(6) = 



fco-l 



k=l 



since we can absorb the correction terms for the slow manifold of the variance into ©(cr^e""). In 
particular, if /cq = 1 then it follows that 



Var(6) = a^Hoiy) + O (a^^e"") . 



(50) 



In principle, we could also calculate higher-order corrections to the slow manifold defined by 
X = H^{Y); see (jl5j, p. 147) and Section[6l For simplicity, we shall only consider the lowest-order 
approximation for a general codimension-two fast subsystem bifurcation. 



For another fast subsystem bifurcation, we will get another approximation of the moments as 
we used ^/—y = x for the slow manifold in the fold scenario. However, we still expect that 

Covins) = AHe{y)] + S2{e,s). (51) 

where H^{y) = YlT=o -^kiv) and 62{e,s) denotes a small e-dependent error term for the second 
moments. In fact, the methods we use here based upon moment equations, integral estimates 
and direct asymptotics all generalize to higher-dimensional phase space and higher-codimension 
bifurcations. Hence we conjecture that (I5ip is still valid for these cases. Although we do not 
calculate the asymptotic relations here, our approach provides a direct computational method for 
the relevant scalings. 

It is very important to recall the result is still only local around the attracting slow manifold 
in a compact set V = V[e). Although (xp,yp) G d'D{0) one always has to use the moment ap- 
proximations by a linearized process outside of a (e, C7(e))-dependent neighbourhood of the critical 
transition point {xp,yp). Small (e, cr(e))-dependent neighbourhoods (R2) near the bifurcation point 
have to be considered separately \17\ [T^. For early-warning signs it is very reasonable to ask 
for the earliest possible statistical indicators. Once a sample path reaches (R2) it is extremely close 
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to a fast jump so a warning sign may be difficult to utilize in applications. 

5 Covariance Scaling Laws near Critical Transitions 

To calculate Ho{y) we have to solve the algebraic equation 

= Ao{y)X + XAoiyf + Fo{y)Fo{yf. (52) 

where the matrices Ao{y) are chosen according to normal form theory from Table [T] (see also 
(j32|) -(|33|) for definitions). It will be convenient to introduce a notation for the symmetric matrix 
Fo(y)Fo(?/)-^ that describes the noise term 

(Ni.iy)) = N{y) := Fo{y)Fo{yf (53) 

for i,j G {1, 2, . . . ,m}. If N{y) is a constant matrix then we deal with purely additive noise while 
dependence on y indicates multiplicative (or slowly parameter-dependent) noise. To distinguish 
between the small noise asymptotics 

e ^ a = a{e) 

and the approach towards the fast subsystem bifurcation point y tending to the origin we use the 
order notation O* for y — )• 0. Recall that (AO) specifies what type of double asymptotics we allow 
and that all results are constrained to a bounded domain i.e. a result w{y) = 0*{W{y)) is to 
be understood as, for a given sufficiently small e > 0, and hence a given o"(e) > 0, there exists a 
compact non-empty domain T)y{e) C with G dT>{Q) but ^^(e) (cf. (Rl) in Figured]) and 
constants ETj, i = 1,2 such that 

KiW{y) < w{y) < K2W{y) 

for all y E T)y{e). In particular, the early-warning signs we are going to derive are for fixed (e, o"(e)) 
sufficiently small, a fixed domain T){e) = T>x{e) x T)y[e) chosen around a slow manifold so that the 
approximation is good as y tends to the transition point but will eventually break down in a small 
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region near the critical transition that scales with e and a and includes the critical transition point 
in its boundary for e = = o". Small (e, (T)-dependent regions containing a critical transition point 
require a special analysis and will not be considered; see the remarks on additional literature in 
Section 8. 

Furthermore, we agree to the convention that any limit as y — t- is always understood as the 
natural one-sided limit if necessary e.g. 0*{^/—y) means y — t- 0^. 

Theorem 5.1. Suppose < e ^ 1 and (A0)-(A4) hold for a fast subsystem bifurcation with one 
fast variable (fold, transcritical, pitchfork, cusp) and e > is sufficiently small. Then the variance 
of the process near an attracting slow manifold approaching the bifurcation satisfies 

Varies) = AHe{y)]+S2{s,e). 

where H^{y) = Ho{y) + 0{e) and 

(VI) (fold) Hoiy) = O; , 

(V2) (transcritical, pitchfork) Ho{y) = O* 

(V3) (cusp) HQ{y) = Oy^ ( ^yt^ ) '" '^^^^^ sZo?« variable y2 multiplies the linear term in the fast 
subsystem normal form ^ . 

In particular, if 62(8, e) <^ a"^ and N{y) is constant then the variance scales, to lowest order, as 
cP' I ^ for the fold, as cP' jy for the transcritical /pitchfork and as cr^/y2 for the cusp transition. 

Proof. We can approximate the variance of the process by its linearization if e is sufficiently 
small. The linearized process has variance 

Var(ei) = a\Ho{y) + 0(e)) + 62(3, e). (54) 

where X = Ho(y) G is the solution of 

= 2Aoiy)X + Niy), ^ x = Hoiy) = (55) 

The non-degeneracy assumptions of the four bifurcations considered are satisfied. By using normal 
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forms, we know from Table [T] that Ao{y) = 0*{^) for the fold transition, Ao{y) = 0*{y) for 
the transcritical and pitchfork transitions while Ao(y) = 0*^{y2) for the cusp transition. Direct 
substitution of these results for ^o(y) into ([55]) gives the result. □ 

The codimension-one fold and the transcritical/pitchfork case in Theorem 15.11 can also be in- 
ferred from previous works see e.g. [I5J. In fact, rigorous proofs without formal truncation are 
available. In these results the higher-order terms do not seem to influence the scaling law in region 
(Rl); this is one of the motivations to consider a formal truncation. The stochastic cusp, and all 
the following codimension-two results, have not been considered previously. It should be noted that 
for fast subsystems with dimension greater than one the stochastic scaling effects are much more 
interesting as the next result shows. 

Theorem 5.2. Suppose < e ^ 1 and (A0)-(A4) hold for a fast subsystem bifurcation with two 
fast variables (Hopf, Bogdanov-Takens, Bautin) and e > is sufficiently small. Then the covariance 
matrix of the process near an attracting slow manifold approaching the bifurcation satisfies 



CoviC,) = a^[H,{y)]+62is,e). 



where H^{y) = Ho{y) + 0{e) and 
(V4) (Hopf Bautin) 



Hoiy) 



( 2 jVii {y)y'^+2Ni2 {y)y+Nii {y)+N22 jy) Nii{y)~N22{y)-2Ni2{y)y \ 

4y{j/2 + l) 4(y2+l) 

Nii{y)-N22iy)-2Ni2{y)y 2N22 (y)y'^-2Ni2 {y)y+Niiiy)+N22 (y) 

\ 4(y^ + l) 42/(2/^+1) / 



In particular, N is a constant matrix with A'^n -|- N22 7^ then 

( 



^o(y) 



V 



(V5) (Bogdanov-Takens; we set 0*{^/—yl) = kyj—yi for AQ{y)) 



Hoiy) 



( -N22{y)+2kNi2{y)^^±2Nii{y)^^+Niik'^yi Nujy) \ 

±4kyi 2 

Niijy) ±2Nii{y)+N22{y)yi/{-yif^'' 

\ 2 2k 
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In particular, if N is a constant matrix and N22 7^ then 



Ho{y) 



O 



\ 



y \yi J 

Nil 
2 



2 

±^ + 0; 



-yi 



Proof. The proof follows the same outline as the proof of Theorem 15.11 Therefore we shall only 
detail the calculations for the proof of (V4). We can again reduce to the linearized process and 
apply the formula 



Cov(e.) = a^Hoiy) + 0(e)) + 62(8, e). 



(56) 



Denote the elements of the scaled covariance matrix as follows 



X = a-^Cov(ei) :- 



Vll Vl2 
Vl2 V22 



Using the normal form matrix ^o(y) from Tabled] we calculate 



Ao{y)X + XAo{yf + N{y) 

-V12 + vuy vi2y - V22 
^ Vll + Vuy V12 + V22y 



+ 



V 



-V12 + Vuy Vll + vi2y 
Vuy - V22 V12 + V22y 



iVii(y) - 2vi2 + 2viiy Nuiy) + vn - V22 + 2vi2y 
y Nuiy) + Vll - V22 + 2vi2y N22{y) + 2vi2 + 2v22y J 



+ 



Niiiy) Ni2iy) 
Nuiy) N22{y) 



(57) 



Equation ([57|) yields three independent conditions. Hence we get a linear system 



2y -2 / Vll 
2y 2 V22 
»^ 1 -1 '^y J \vi2 J 



^ I -Niiiy) ^ 

-N22{y) 
V -Nuiy) J 



that can be solved for (vii, U22, ^^12)- This result can be substituted a.s X = HQ{y) in ([56]) and this 
yields the first part of ( V4) . If is a constant matrix then direct asymptotics shows that 



Nii+N22 _^,fl 

4y y \y 



N11 + N22 ,^^fl 
4^ = [y 
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where y — J- as we approach the critical transition via the attracting slow manifold. For the 
covariance we get 



Nu - N22 2Ni2 Nil - N22 , , , 
vi2 ~ —y = + Oy {y) 



The result for the Bogdanov-Takens transition follows by the same techniques. □ 

Before we continue to codimension two bifurcations in M? and M^, let us interpret the results of 
Theorem 15.21 for 52{s^ e) <^ o"^. For the Hopf transition with a fixed noise level ex > we have found 
that the variance of the coordinates increases as 0*{\/y) as the bifurcation point is approach with 
y — )• 0. This result is expected as we already saw an increase in variance for the one-dimensional 
fast subsystem bifurcations. However, for the covariance the additive noise case with a constant 
matrix N yields 

This implies that the covariance tends to a constant as y — t- 0; even more surprisingly, for the 
reasonable assumption of equal individual diffusion Nu = N22 we get that the covariance tends to 
zero as the bifurcation is approached. Hence we can already conclude that measuring covariances 
can also provide important information to predict critical transitions. For the Hopf transition with 
multiplicative noise, let us just consider the simplest case of linear multiplicative noise without 
correlation Nii{y) = ciy, N22{y) = C2y, Ni2{y) = 0. Then we find 

Var(6,s) ^^^^-^-=^ = 0;(l)+0;(/), ifci/-C2, 



(ci + C2)y 


2ciy3 


Ay 


y 


(ci + C2)y 


2c2y3 


Ay 


y 


(ci - C2)y 


= Ol{y) 


4 



Var(6,s) ^^:^-^-^ = 0!(l) + 0!(y2), ifci/-C2, 



Cov(^i,s,6,s) 



Therefore measuring the variance alone is not expected to yield valuable information; indeed, 
variance tending to a constant could be interpreted as a normally hyperbolic regime without critical 
transitions for additive noise [65] . Obviously one could discuss further interesting scalings depending 
on the matrix N[y). It should be clear from the formulas (V1)-(V5) and the previous discussion 
how to approach these situations as long as the system is in normal form near the bifurcation point. 
We proceed to look at some results for the remaining codimension two bifurcations. 
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Theorem 5.3. Suppose < e ^ 1 and (A0)-(A4) hold for a codimension-two fast subsystem 
bifurcation with at least three fast variables (Gavrilov-Guckenheimer, Hopf-Hopf). Assume that e 
is sufficiently small and that N = N{y) is a constant matrix. Then the covariance matrix of the 
process near an attracting slow manifold approaching the bifurcation satisfies 

Covins) = AHeiy)]+S2{s,e). 



where H^{y) = Ho^y) + 0{e) and 

(V6) (Gavrilov-Guckenheimer) if Nn ^ and N22 + -^33 7^ then 



Ho{y) = o* 



( _j_ 

Ni3 



Ni3 

UJ 

y2 



\ 



N12 N22-N33 



N12 

N22-N33 
J_ 

y2 



(VI) (Hopf-Hopf special case: Ao{y) given by i^) if Nn + A^22 / and N3S + iV44 / then 



Ho{y) = 01 



\ 





NX1-N22 


Ni4,UJ2-N2-il^l 


Af24'^l-A''l3'^2 


y\ 






i i 


N11-N22 


J_ 


NiiUj-i+N2ii^2 


NuLUi — N2-iOJ2 




y\ 


? ^ 




Ni4,U2~N23Ull 


Nri.'jJi+N24,'^2 


J_ 


7V33-7V44 


^i-^2 




y2 


4aJ2 


N2AUJ1-N-13UJ2 


Ni4Ldi — N23t02 


N33-NAA 






ujI-uI 


4£J2 


2/2 



We shall omit the calculations for the proof of Theorem 15.31 as it follows the same steps as the 
proofs of Theorems 15. 1115. 2[ We remark that the solution of the algebraic equation (j52p becomes 
much more cumbersome for systems in M'^ and and we compared our solution to the results 
obtained by a computer algebra system [53j. Another important note on Theorem 15.31 is that we 
do not have to assume explicitly that uji 7^ uj2 for the Hopf-Hopf bifurcation since this is included 
in assumption (A2). The 1:1 resonance case at a Hopf-Hopf bifurcation wi = uj2 (see \9A\ 143] ) 
naturally appears as a special case in our analysis. In particular, when \ujf — a;|| is small then 
the covariances of the two-by-two off-diagonal blocks in (V7) can also get large near a Hopf-Hopf 
critical transition. 
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6 Double- Singular Variance Asymptotics for the Fold 

In the last section we have computed the leading-order term for the covariance near a critical 
transition for all fast subsystem bifurcations up to codimension two. In current applications of 
critical transitions one frequently encounters the fold bifurcation. Prom a mathematical viewpoint, 
the fold bifurcation is the lowest slow codimension, lowest fast dimension generic bifurcation without 
further assumptions. Both reasons warrant a more detailed asymptotic study to determine higher- 
order correction terms to the formula 



Var(e.) = 



+ 



from Theorem 15.11 for additive noise. We can always assume a preliminary normal form transfor- 
mation |731 [69] and consider on the slow time scale 

^ (58) 

dy = —1 ds, 

where we assume additive noise to simplify the algebraic manipulations to follow. We also refer 
to Figure [U and consider the attracting branch Cg = {{x^y) € : x = y/y = h^iy)} of the 
critical manifold. Fenichel's Theorem provides an attracting slow manifold = {{x,y) G : x = 
h^{y) = hQ{y) + 0{e)}. A converging series expansion for Cq can easily be derived by direct regular 
asymptotics where convergence of the asymptotic series is guaranteed by Fenichel's Theorem. 

Lemma 6.1. The attracting slow manifold for the fold bifurcation normal form is given by 

Terms of order 0{e^) or higher are omitted but can easily be calculated from a recursive solution of 
algebraic equations. 

Setting = Xs — h^iy) and using Lemma 14.21 for equation (I58p we find that the scaled variance 
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Xs = cr-2Var(^s) satisfies the ODE 



eX = -4:hJy)X + l, 

^ ' (60) 

y = -1. 

The attracting critical manifold of ([60]) is given by Co = {{X,y) £ M."^ : X = Ho{y)}. Fenichel's 
Theorem yields an associated attracting slow manifold = {{X,y) £ M? : X = H^{y) = HQ{y) + 
0{e)}. We already know from the proof of Theorem 15. II that Ho^y) = l/(4y^). 

Proposition 6.2. The attracting slow manifold associated to (f60ll has an asymptotic expansion 
given by 

^ 1 3 n 7 , 201 4 3837 , , 

= + ^3V + + ^ W + ^ + ^''^ 

Terms of order 0{e^) or higher are omitted but can easily be calculated from a recursive solution of 
algebraic equations. 

Proof. We make the ansatz H^{y) = H^iy) + eHi{y) + e^H2{y) + • • • . Using this ansatz and the 
result from Lemma l6.1l in (I60p we get a hierarchy of algebraic equations at different orders 

= I - Aho{y)Ho{y) 
^ = -4 ^ H.iy)hM 

i+j=k 

where k £ {1,2, . . .}. The result (I6ip follows by direct calculation. □ 

We expect that the expansion up to fourth order of is sufficient for all practical purposes. 
Recall from the end of Section [4] that the condition e'^°'^"i7fco(y) = 0{a) determines whether terms 
of the expansion for can be moved to the higher-order correction 0{a^/e°'). Proposition 16.21 
yields the conditions 

0{a) for /CO G {1,2,...}. (62) 



y(3fco+l)/2 

For the fold bifurcation, we know that the critical scaling of y to stay inside the normally hyperbolic 
regime is y ~ e^/^; see Section [2] and assumption (AO) as well as Lemma |6. 11 Suppose y ~ e^" for 
some a < 1/3 and use the scaling in (j62|) for A;o = 1 then e^^^a _ 0(u) ig the condition to move 
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all slow manifold correction terms into the higher-order for the variance estimate; for a = this 
condition obviously reduces to the known fact e = 0{a) from equation (f50]) . 

Using the critical scaling 2a = 2/3 in ()62p we get the condition e^~^ = 1 = 0{a) which can 
never hold under assumption that a{e) — )• as e — )• by assumption (A4). Therefore, the slow 
manifold approximation in (Rl) obtained by the linearized process for the moments is not valid 
in (R2). 

7 Applications 

We are going to present five applications to illustrate the previous results. We also indicate how 
novel conclusions about the applications follow from the theory. 

7.1 A Climate Box-Circulation Model 

The Stommel model [88] describes the North Atlantic Thermohaline Circulation (THC) by 

two boxes Bi and B2 representing low and high latitudes respectively. An atmospheric freshwater 
flux and differences in insolation can induce temperature and salinity differences AT = T1 — T2 and 
AS" = Si — 82- The resulting system has an Atlantic northward surface current and an Atlantic 
southward bottom current. For a version of Stommel's box model |24] it can be shown that 
(AT, AS") obey a two-dimensional fast-slow system where the temperature difference represents the 
fast variable pjj. After reduction to an attracting slow manifold and a re-scaling of the variables 
the dynamics reduces to 



where Y represents the salinity difference, we fix r/^ = 7.5 and /i is a parameter proportional to 
the atmospheric freshwater flux. Obviously the freshwater flux can also be viewed as a dynamical 
variable and we assume that it changes slower than Y. Furthermore we assume that (|63p is subject 
to small stochastic perturbations which is reasonable if we decide not to model the system in more 
detail. Setting x := y and y := fi we get another two-dimensional fast-slow system 




(63) 



dxs 



i [Vs - Xs{l + 7.5(1 - Xs f)] ds + ^^F{ys)dW, 



(64) 



dys 



g{xs,ys)ds. 
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The deterministic critical manifold is Cq = {(x, y) £M? : y = x(l + 7.5(1 — x)^) =: ho^x)}, which is 
immediately recognized as a classical S-shaped (or cubic) fast subsystem nonlinearity. There are 
two fold points (fast subsystem fold bifurcations) at 

(.-,,-)= (-L(10-v/l5),ll + -i=) a„d (1(10 + v/l5),il--i=)^ 

The critical manifold splits into three parts Cq' := Cq Ci {x < x-}, Cq := Co H {x~ < x < x~^}, 
and Cq'~^ := Co n {x > x^"} where Cq'^ are attracting and Cq is repelling. The lower branch Cq'~ 
represents small salinity difference which corresponds to a weak THC. The upper branch Cq'^ 
corresponds to a strong THC which can be viewed as the present state of the climate. A critical 
transition from a strong to a weak THC would mean a significant cooling of the mild European 
climate. Therefore, we shall focus on the critical transition (x"*", y"*") with initial conditions on Cq'^. 
The initial condition will be fixed at (xo,yo) = ixo,3/2) £ Cq'"*" which roughly corresponds to the 
drop point [60J on the upper attracting critical manifold after a transition at (x^,y^). 
We start by simulating (j64p using an Euler-Maruyama method [49j using 

e = 0.01, a = 0.01, F{y) = l, g{x,y) = -I. (65) 

where the assumptions on g mean that one may also interpret y as a time variable. A typical 
sample path is shown in Figure [IJa); the path is stopped at a final value y = 0.95. We want to 
estimate the variance Var(ys) from a time series 

yo = yso'^si, • • • ,ysiv = 0-95, x^^jX^,, . . . ,Xs^. (66) 

The values be viewed as functions of y since = (s — sq) + yo and we indicate this by 

writing Var(x(y)) := Var(xs). The goal is to estimate the variance. There are several possibilities 
to extract an approximation: 

(Ml) Consider a single time series. Select a moving window of fixed length M and compute the 
sample variance for M + 1 consecutive points xj, . . . ,Xj-^M] see Figure [2]^a)-(b). This pro- 
vides an estimate for the variance Var(x(y)) roughly at the midpoint of the moving window 
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1 .025 ' — 
1.04 
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1.02 



50 



100 



Figure 2: Illustration of the different techniques (M1)-(M4) to approximate the variance Var{x{y)); 
we use the Stommel-Cessi model (j64p with parameters given in (j65p . (a) Typical time series (black) 
near the attracting critical manifold C"'"*" (red) up to the fold point {x~^,y~^) (black dot). We also 
show two sliding windows (green) where the dashed green line is a linear trend and the solid green 
line is given by Cg'^ i.e. the green curves are used for linear and CM detrending respectively, (b) 
Detrended time series xd from (a) corresponding to the two (green) sliding windows, (c) Zoom 
near y = 1.05, 5 sample paths are shown. The dots (magenta) mark the five points of the paths at 
y = 1.05. To calculate the variance Var(x(j/ = 1.05)) one simulates many paths, (d) Simulation of 
the fast subsystem of (j6l|) with y=1.05 for a fixed fast time t G [0, 100]. 

Jd Sfc^o Vj+k- The idea is that if the window is sufficiently small and we have sufficiently many 
data points inside each window then we can calculate a good approximation to Var(x(y)) for 
each y. 



(M2) Consider a single time series as for (Ml). We can remove a given trend from (|66p before 
calculating the variance. For example, interpolating (j66p linearly and subtracting the resulting 
linear function from the time series yields a variance estimate with linear detrending. 
Another natural possibility is to remove the critical manifold as a trend; we call this critical 
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manifold (CM) detrending. See also Figure [2]^a)-(b). 

(M3) Another possibility is to consider a large number R of time for 
r G {1,2, and then calculate the variance Var(x(yj)) at yj as the sample variance 

of {x^^\x^^\ . . . , x^^^}. This idea is illustrated in Figure [2]^c) and avoids the moving window 
technique. However, it does require multiple time series passing near the same critical point. 

(M4) Instead of simulating the entire SDE ()64p we can also assume that y = yj is constant, simulate 
the fast subsystem for a sufficiently long time and then calculate Var(x(yj)) from this fast 
subsystem time series; see Figure [2fd). 

Each of the methods (M1)-(M4) has different advantages and disadvantages. A direct sample 
variance measurement using the sliding window technique (Ml) does include the curvature of the 
critical manifold in the estimate as demonstrated in [65]. Linear detrending requires no a priori 
knowledge about the dynamics but can obviously not remove curvature near the fold point. CM 
detrending corresponds to the change of variable '■= Xg — hQ{ys) which is closest to the theoretical 
situation discussed in Sections[l][6l However, this requires a priori knowledge of the critical manifold. 
The method (M3) requires many sample paths which is a restriction while the method (M4) requires 
the ability to simulate/measure the fast subsystem for a long time. In Figure [3] we compare the 
different methods for the Stommel-Cessi model. Figures [3|^a)-(e) provide the variance estimates 
together with a least squares fit of 

Var(x(y)) = jj_ (67) 

with fitting parameters A and yc- The results in Figure [3] show that all methods can capture 
the variance increase as predicted by the theory. The sliding window technique seems to deviate 
the most from the theory compared to the other four methods but it requires the least amount 
of data as one basically produces a plot similar to Figure [3^a) with just a single time series. By 
fitting ()67p we also obtain an estimate for the critical transition point yc which is slightly delayed 
due to positive e. All techniques capture this effect. We get the estimate that yc G [0.92, 0.95] 
which is a very good prediction compared to direct simulations. Overall one may conclude that the 
theoretical predictions of variance increase near a fold point apply very well in the context of the 
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Figure 3: Comparison of different methods to estimate the variance V = Var(x(y)) for the Stommel- 
Cessi model (j64p with parameters given in (j65p . The black curves in (a)-(e) indicate the variance 
estimate and the green curves are obtained by least squares fit of (I67p . (a) Sliding window technique 
(Ml) without detrending, average over 1000 sample paths, (b) Sample paths "pointwise variance" 
(M3), average over 1000 sample paths, (c) Sliding window with linear detrending (M2), average 
over 100 sample paths, (d) Sliding window with CM detrending, average over 100 sample paths, 
(e) Fast subsystem simulation (M4) for a fast time t G [0, 100]. (f) The critical manifold (red/blue) 
with fold point (black) is shown. The green markers indicate the estimators for yc from a least 
squares fit of (j67p plotted at the same x- value as the fold point; the green star "*" is the lower 
bound estimate for yc from (a) and (e), the green circle "o" marks yc for (b), the green plus "+" 
corresponds to (c) and the green "x" marks yc for (d). 



Stommel-Cessi model (j64p and that the different time series analysis methods all have advantages 
as well as disadvantages depending on the situation. Obviously we do not make any claims about 
the real THC with our calculations as this requires the analysis of temperature data sets. 



7.2 Epidemics on Complex Adaptive Networks 

Consider a network of social contacts and a disease that can be spread via these contacts as 
described in [H]. Individuals of the population correspond to nodes (or vertices) and social 
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contacts correspond to undirected links (or edges). Denote the total number of nodes by 
and the number of hnks by K and assume that and K are constant; define the mean degree 
^ := 2K/N . The dynamical states of the nodes are either susceptible or infected giving classical 
SIS dynamics. If a link between an infected and a susceptible node exists then the susceptible 
node becomes infected with probability p at each time step. Infected nodes recover to susceptible 
status with probability r. Susceptibles might try to change their connection from an infected node 
to a susceptible one. To model this effect we assume that the network is adaptive so that the 
topology of the network influences the dynamics of the nodes and vice versa. We use the following 
dynamical variables to describe the SIS adaptive network 

Xl := # {mfcctcd} _ "(jgnsity of infected individuals" , 

# {links between infected and infected} -x j -j. r tt i- i j) 

X2 := jY = per capita density of ii-lmks , 

# {links between susceptible and susceptible} u -x j -j. r oo i- i » 

X3 := = per capita density of Sb-lmks . 

Note that the density of susceptible individuals is (1 — xi) and the per capita density of Sl-links 
is (/i/2 — X2 — X3). To capture the full adaptive network dynamics one would have to take into 
account also triples (triangle subgraphs) and all other higher-order (network) moments. We use 
the moment closure pair approximation [59] to express the higher-order moments in terms of x 
which yields [H] 

x'l = - X2 - X3) - rxi, 



x 



= p{!±-x2-x,)(^2^^ + lj-2rx2, (68) 
x'3 = [r + w){l^ - X2 - X3) ^ 

For our analysis we fix the following parameters 

r = 0.002, w = 0.4, N = 10^ = 10^ ^ = 20. (69) 

Assume that p is a slow variable and increases over time. For example, we could think of a virus that 
evolves towards a more infectious variant in time. Using the standard notation for slow variables 
we let y := p and assume y' = e. It is also reasonable to consider the scenario that the density of 
infected nodes and the link densities can exhibit stochastic fluctuations; in particular, this might 
lead to a model that is more realistic than the moment closure ODEs. Combining this assumption, 



41 



the slow equation and ([68]) we get 



dxi 

(1X2 
dX3 

dy 



X2 
X2 



(r + u')(f 



X3) 

X3) 

- X2 



rxi] ds + ^^dW'^'^\ 



-X2-X3 



l—xi 



X3) 



+ 1 



l — Xl 



ds + ^JW^'^\ 
ds + ^dW^^\ 



(70) 



1 ds, 



where we omit the subscript s for Xg and Wg = {Ws^\wP\Ws'^^)'^ for notational convenience. 
Although the algebraic expression for the deterministic critical manifold Co of (j70l) can be computed 
we shall only focus on the subset 



Co* := [{x,y) G ([0, 1] x [0,///2]2) x [0, 1] : = = rro, X3 = |} C Cq. 



The solution xi = = X2 and 2:3 = ij,/2 corresponds to an equilibrium point of (|68[) with no infected 
nodes that can also be obtained by considering the initialization of the network as a random graph 
[41j . The fast subsystem linearization around Cg is given by 



Dxf\ 



-r -y -y ' 

—y — 2r —y 
y y/j, — r — w y^ — r — w J 



(71) 



Using the parameter values (j69p and (j7ip we can easily calculate that a single eigenvalue of 
(j71|) crosses the imaginary axis at y = yc = 0.0201. Another direct calculation shows that Cq splits 
into two subsets Cq" = {y < yc} H Cq and Cg'" = {y > yc} H Cq where Cg" is normally hyperbolic 
attracting and Cq*" is normally hyperbolic repelling. Note that the fast subsystem bifurcation of 
the trivial solution Cq to (|68|) suggests a transcritical or a pitchfork bifurcation. In Figure S] we 
show part of the critical manifold Cq including the trivial solution Cq; the computation has been 
carried out using numerical continuation [39j . Figure H] shows that the bifurcation is transcritical 
and y = yc is the infection probability threshold. 

For direct simulation of ()69p one has to ensure that x G [0,1] x [0, /u/2]^ as the densities 
are constrained. Therefore, we set a point that lands outside of the domain at a time step to its 
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0.02 0.04 y 0.06 0.02 0.04 y 0.06 0.02 0.04 y 0.06 

Figure 4: Parts of the critical manifold Cq for the SIS-model (I70p where attracting branches are 
red and repelling branches are blue; parameters are given by (|69p . The manifolds (fast subsystem 
equilibrium branches) have been computed using numerical continuation [39]. A transcritical bi- 
furcation (branch point, [BP]) is detected at y = yc = 0.0201. For the number of infected nodes 
we show the continuation of Cq away from the branch point; it undergoes a fold bifurcation (limit 
point, [LP]) and stabilizes at a supercritical Hopf bifurcation [Hj. 

associated boundary value, e.g. if xi{sj) < for some numerical time step Sj then we set xi(sj) = 0. 
This simulation is formally outside of the theory developed in Sections [2][6l Nevertheless, Figure [5] 
shows that the theoretical results are useful. Figure [5]^a)-(c) shows a typical sample path and we 
see that the xa-coordinate in (c) starts to decrease beyond the singular limit critical point whereas 
the other two variables do not show any recognizable trend in (a)-(b). It is interesting to note that 
the density of infected individuals does not seem to play a role as an early-warning sign for the 
epidemic outbreak. 

Figure [5)^d) shows the variance V = (^4,^2,^3) = (Var(xi), Var(x2), Var(x3)) associated to 
the sample path in (a)-(c) by using a sliding window technique; the size of the sliding window 
corresponds to the gap in the curves near y = 0. Figure ^e) shows an average variance Vi for 
i = {1,2,3} over 1000 sample paths. Observe that X3 is the best predictor variable and this leads 
to the conjecture that the increase in variance should scale like the inverse of the distance to the 
critical transition; see Figure [5]^f). Note that the critical transition aX y = yc for e = is delayed 
due to the time scale separation [65] . We conclude from our results that it is crucial what property 
of a complex system we actually measure to make predictions. Indeed, the SIS-epidemic model 
suggests that measuring the variance in links can be much more important than just the number of 
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Figure 5: Simulation results for ()70p with boundary conditions to constrain x = (xi,X2,a;3) € 
[0, 1] X [0, ///2]^; parameter values are given in ([69]) and (cri, cr2, cts, e) = (0.01, 0.01, 0.01, 0.005). (a)- 
(c) show a time series and (d) shows the associated variance of this series calculated by a sliding 
window technique, (e) Average V of the sliding-window variance for 1000 sample paths; we see that 
V3 shows an increase near the bifurcation, (f) Inverse of averaged variance \/V where we clearly 
see that V3 scales like l/(yc ~ v) up to a delayed epidemic threshold y%. We also show two linear 
fits to (Va)"^, one before the threshold (early-warning regime) and one after the threshold (start of 
critical transition). The actual full epidemic outbreak is not shown in the plot and occurs roughly 
between y = 0.05 and y = 0.07. 

infected individuals. Furthermore, the technique we developed here can also be applied to adaptive 
networks in completely different contexts [42^ [67] . 

7.3 A Switch in Systems Biology 

To understand complex molecular networks one often seeks to construct models of simpler building 
blocks of the network. These building blocks are composed of genes and proteins and can often act 
as various kinds of "switches" inside a more complex system. Dynamical systems methods for these 
systems biology questions are a highly active research area [20] . Low-dimensional dynamical systems 
have been proposed to model the smallest units in a molecular network. A typical example is the 
activator-inhibitor system. Suppose the activator species R is produced in an autocatalytic 
reaction but rising R also promotes the production of an inhibitor species X. More concretely, 
one may think of both species {R, X) as concentrations of proteins. Activator-inhibitor systems 
incorporate positive and negative feedback which can lead to oscillations. One model proposed for 
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activator-inhibitor oscillators [93] is 



R' = koG{k3R,k4,Ji,J2) + kiS-k2R-k7XR ^^^^ 
X' = k^R-k&X 



where the Goldbeter-Koshland function G [38^ 178] is 

2uK 



G{u,v,J,K) 



V — u + vJ + uK + y^{v — u + vJ + uK)'^ — 4(u — u)uK 



and kj for j G {1,2,3,4,5,6,7}, Jj for i G {1,2} and S are parameters. The main bifurcation 
parameter is the signal strength S which can be viewed as an external input to the system (j72p . 
We are going to fix the other parameters following |93j as 

kQ = 4, ki = k2 = k3 = k4 = kr = 1, k^ = 0.1, ke = 0.075, = J2 = 0.3. 




Figure 6: Dynamics for e = for the deterministic version of the activator- inhibitor system (173p . 
The critical manifold Cq is the red-blue curve which looses normal hyperbolicity at a fast subsystem 
subcritical Hopf bifurcation (black dot, [H]) at y ^ 0.09146. The generated small limit cycles (blue) 
are first repelling and then undergo a fold (or saddle-node, or limit point [LPC]) bifurcation; the 
large fast subsystem limit cycles (red) are attracting. A critical transition occurs near the Hopf 
bifurcation as trajectories leave the critical manifold and jump to a large limit cycle. See also 
Figure [7| for the fast subsystem phase portraits. 
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Xl Xl Xl 

Figure 7: Illustration of the subcritical Hopf bifurcation for the deterministic fast subsystem of 
(f73l) : equivalently the results apply to ([72]) with [R, X) = (xi, X2). Nullclines are shown in magenta 
for Xl and in orange for X2- Trajectories are black and invariant sets are red (stable) and blue 
(unstable), (a) y = 0.01: The system has a stable spiral sink, (b) y = 0.085: In addition to the 
spiral sink there exist a small unstable and large stable limit cycle, (c) y = 0.12: The equilibrium 
point is a spiral source and only the large stable limit cycle exists. 

Let us consider the case when the external input S" is a slow signal that starts out sufficiently 
low so that no oscillations occur for (I72p . Then we let S =: y increase until a transition to 
large oscillations is observed. It is reasonable to assume that the variables {R,X) =: {xi,X2) are 
stochastic with correlated noise. Under these assumptions we can write (I72p as the SDE 



dxi = i[4G(xi,l,0.3,0.3)+y-xi-xiX2]ds + -^(FiidII^(i)+Fi2dVF(2)), 

dx2 = i[0.1xi- 0.075x2] ds + ^(F2idiy(i)+F22(iiy(2))^ (73) 

dy = 1 ds. 

The critical manifold Co for the deterministic part of ()73p is given by 

Co = !^{xi,X2,y) G : X2 = ^2:1, y = xi + ^x? - 4G(xi, 1, 0.3, 0.3) 

It is easy to check that the critical manifold is attracting for y < yn,! and y > yH,2 and repelling for 
yH,i < y < yH,2 where yH,i ~ 0.091462 and yH,2 ~ 0.440903 are fast subsystem Hopf bifurcation 
points. We focus on the subcritical Hopf bifurcation at y = yn,!- Figure [6] shows an illustration of 
the singular limit dynamics near this Hopf bifurcation point. Repelling fast subsystem limit cycles 
are generated at the Hopf bifurcation. These cycles undergo a further fold (or saddle-node, or limit 
point) bifurcation to attracting cycles which grow rapidly. By looking at the phase plane of the 
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fast subsystem in Figure [7| we observe that the xi-nuhchne can also be viewed as another critical 
manifold of the two-dimensional system {xi,X2) where X2 would be fast and xi be even faster which 
yields a three-scale system with canard explosion [61j. The important outcome of this mechanism 
is that passing from the attracting critical manifold 

:= {{xuX2,y) GR^:y< y^.i} n Co 

through the Hopf bifurcation produces a critical transition to large limit cycle oscillations. The 
critical transition can be viewed as an almost instantaneous switch to sustained oscillations. For 
the stochastic simulation we recall from the definition in equation (j53p that 
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For numerical simulations fix A'^n = 1 = N22 and A'^12 = 0.2. In Figure [8|^al)-(a3) the variance 
and covariance near the subcritical Hopf bifurcation aX y = yn,! for the activator-inhibitor system 
([73l) are shown. The variances Var(xi^2) have been fitted using 

Var(xj(y)) = , for j G {1, 2} (74) 

y Uc 

with fit parameters A and yc- The covariance has been fitted linearly. The variance of the fastest 
variable \ax{xi{y)) behaves approximately as predicted near the critical transition as 0{\/y). 
However, the variance Var(x2(y)) of the slower variable X2 does not show a clear increase and 
the covariance near the critical transition is not constant. This shows that the three-time scale 
structure requires a very careful analysis and a transformation to normal form would be needed to 
apply Theorem 15. 21 A prediction of the critical transition point can still work e.g. using Var(xi(y)) 
produces the estimate yc ~ 0.094. We have also compared the activator-inhibitor results to a Hopf 
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Figure 8: The row labels denote Vi = Var(xi), V2 = Var(a;2) and Ci^2 = Cov(xi,X2). (al)-(a3) 
Parameter values are e = 10"^ and a = 10"^ for ([73]). (bl)-(b3) Parameter values are e = 5 x 10 
and a = 10^^ for (j75|) . All figures have been computed from 100 sample paths by a sliding window 
technique (black curves). The variances have been fitted using (j74p and the covariance have been 
fitted linearly (green curves). We observe that the normal form corresponds perfectly to the theory 
but that the three-time scale structure of the activator-inhibitor system becomes visible in the 
variance and covariance measurements. 

bifurcation normal form system 

dxi = I [yxi -X2 + XI (xl + xD] ds + ^^ {FudW^^'^ + FisdW^^^)) ^ 

dx2 = i [xi + ya;2 + X2(x2 + a;2)]ds + ^(F2idl^« + F22dVF(2)), (75) 

dy = 1 ds. 

Figure [8)^bl)-(b3) shows the results which match Theorem 15.21 as expected. For the covariance 
there is a clear difference between (|75p and the activator-inhibitor system; compare Figures [8)^a3) 
and[8]^b3). The increase of the covariance near the critical transition is not expected and might be 
related to deterministic rotation around the slow manifold i.e. the manifold is attracting but 
also a spiral sink of the fast subsystem; see also Section 17.41 where another possible explanation is 
given. 
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To conclude, observe that the bifurcation structure displayed by (j73p has a fast subsystem with 
an S-shaped (cubic) critical manifold which makes the results applicable also to typical neuroscience 
models such as bursting neurons [54:\ I82j . Therefore, we have shown that subunits of molecular 
networks as well as neurons in neural networks do have information available that allows them to 
predict a future state without previous knowledge of the exact position of this state. Whether this 
predictive potential is actually used in a real molecular or neural network is far beyond the scope 
of this paper but certainly constitutes a fascinating question. For a recent application to excitable 
neuron models and epileptic seizures see [72] . 

7.4 A Predator- Prey Systems near Codimension Two Bifurcation 

Sudden shifts in ecosystems have been a primary motivation to develop the theory of critical 
transitions |85) . Recently also experimental evidence has been provided [33]. However, many 
studies seem to view fold critical transitions as the only relevant transition [US]. This viewpoint 
does not seem to be appropriate as codimension two (and higher codimension) bifurcations occur 
very frequently in ecological models [8]. Here we focus on the analysis of a classical predator-prey 
model [8] 

where x\ represents prey, X2 represents predators and a, (5, ^, 7 are positive parameters. The 
bifurcation analysis of (|76|) in the (a, 5)-parameter plane has been nicely described by Kuznetsov 
(see [69j, p. 327-332) under the assumptions 7 = 1 and < ^ ^ 1. For numerical simulation we fix 
7 = 1 and ^ = 0.01. 

We set y\ := a and 7/2 '■= S to indicate that these parameters will be viewed as slow variables. 
Part of the bifurcation diagram for ()76p is shown in Figure [9l Figure [9] shows two curves of fold 
bifurcations, which actually form a closed curve clp ("isola") in parameter space. This curve has 
a tangency with a supercritical Hopf bifurcation curve ch at a codimension-two Bogdanov-Takens 
(BT) point. We do not show the homoclinic bifurcation curve originating at the BT point in Figure 
[9l The curves clp and ch can be calculated explicitly [69]. One simply uses the linearization of 
DxF{x*) of ([75]) at an equilibrium point (xi,X2) = x* and applies the conditions det^D^F^x*)) = 
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Figure 9: Partial bifurcation diagram of the Bazykin predator-prey model ()76p with 7 = 1 and 
^ = 0.01. The parameters (2/1,2/2) can be viewed as slow variables. The main organizing center 
in the diagram is the co dimension- two Bogdanov-Takens (black dot, [BT]) point that occurs at a 
tangency of Hopf (red, [H]) and fold (blue, [LP]) bifurcation curves. Phase space diagrams for the 
different regions Qi, Q2 and Q3 are shown in Figure [TOl note that splits into two sub-regions 
by a homoclinic bifurcation curve which we do not show here. The dashed curve (green) shows a 
slow subsystem trajectory that approaches the BT point. 




4 8 12 4 8 12 048 12 



Xl Xl Xl 

Figure 10: Phase space diagrams for different parameter regions in Figure [U black curves are 
trajectories. Qi: (2/1,2/2) = (0.45,0.35); Q2: (2/1,^2) = (0.35,0.3); Q3: (2/1,2/2) = (0.45,0.15). In 
Qi there is a stable spiral sink outside of the chosen range at {xi,X2) ~ (92.12,3.34). A spiral 
sink equilibrium point also exists in Q2 and outside of the displayed ranges. In Q2 we have a 
spiral sink (red dot) and a saddle point (blue dot) that correspond to attracting and saddle-type 
branches of the critical manifold. In Q3 we have a spiral source and a saddle point corresponding 
to unstable and saddle- type critical manifolds. 

and Tt{DxF{x*)) = 0; this gives 

CLP = {2/GK2: 4e(yi _ 1)3 + ((y2 _ 20y, _ 8)^2 + 22/1^(2/2 _ llyi + 10) 

+2/f (2/1 -l)')2/2- 4(2/1+^3^2^0}, 
CH = {y£R^: 4^(2/1(2/1 - 1) + e(2/i + 1)) + (2(C + 1)2/? + (3^'- 2^-1)2/1 

+ae - 2^ + 5))2/fo+ (2/1 + e - i?yi}. 



The Bogdanov-Takens point satisfies all genericity conditions required by assumption (A2) so that 
Lemma 13.71 and Theorem 15.21 apply. The normal form coefficient is s = — 1 in equation (jl6p . Hence 
the only stable equilibrium point near the BT-point can be found between the Hopf and fold curves 
in region Q2 in Figure [9l Phase portraits for different regions are shown in Figure [TOl A stable limit 
cycle can occur between the Hopf and homoclinic bifurcation curves in region but we ignore this 
possibility and restrict ourselves to critical transitions via fast subsystem stable equilibrium points. 
It is natural to assume that the populations (xi,X2) are subject to stochastic fluctuations and 
to view (2/1,1/2) as slow dynamic variables, changing slowly due to evolutionary or environmental 
effects. This converts (1761) into the SDE 



dxi = - 



_ XIX2 _ p Z 
-^1 1+axi ^-^1 



dyi = gi{x,y)ds 
dy2 = g2{x,y)ds 



(77) 



where we have assumed uncorrelated noise in the fast variables. The critical manifold Co of the 
deterministic version of (|77p has an attracting branch Cq in the region Q2 (see Figures l9l and ITOj) 
corresponding to a spiral sink of the fast subsystem (f76|) . We want to approach the Bogdanov- 
Takens critical transition via a slow flow inside the region Qz- Figure HUl shows a dashed curve 
(green) which is a possible slow subsystem trajectory. It is part of a candidate 70 that undergoes 
a critical transition according to Lemma l3.7[ In principle, we could try to embed such a candidate 
into an explicit slow flow y = g{x,y) = {gi{x,y), g2{x,y))'^ . 

For numerical simulations of (|77p it will suffice to define a single trajectory 70 along which 
we approach the BT transition. We can obtain 70, for example, by polynomial interpolation 
of a suitable set points lying in Q2 and the BT point. The initial condition for our numerical 
simulation is chosen as (xi, X2, yi, 2/2) ~ (3.1544, 1.8849,0.3,0.3293), where the y-coordinates lie on 
the dashed curve indicated in Figure El and the x-coordinates are on the attracting critical manifold 
Cq. Calculations have been carried out for 50 sample paths and the variance has been calculated 
via a moving window method for each path (see the gap in Figure ITlTa) for the window size) with 
linear detrending. Then the results is averaged over the 50 paths. Figure ITTT a) compares the 
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variances Vi = Var(xi(y)) for i E {1,2}. 



Fits of 
Vi = Var(xi) i ^1 



X 10' 
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Figure 11: Simulations averaged over 50 sample paths of the Bazykin predator-prey model (|77p 
with 7 = 1,^ = 0.01 and (e, a) = (3 x 10"^, 1 x 10~^). (a) Variance curves Vi = Var{xi{y)) for 
i G 1,2; the red curve corresponds to Vi and the black curve to V2. In (b) and (c) we repeat these 
curves and show different fits. The green curves correspond to (fTHj) and the blue curve to ([79]) . 



Figures [TTTb)-(c) show fits (green curves) of the variances 



Var(xi(y)) 



A 



yi,c - yi 



for i G {1,2} 



(78) 



and also an inverse square-root fit (blue curve) 



Var(x2(y)) 



A 



\fy^,- 



2/1 



(79) 



where A, yi^c are the fitting parameters . Note that both variances increase like O* (l/(yi,c— yi)) near 
the critical transition and that (j78p is a good fit for V2 while (|79p is not. At first, this might look 
unexpected since the normal form analysis predicts one variance to increase like Cy(l/\/yi,c — yi)- 
However, equation (j77p is not in normal form. To explain the effect let us consider the Bogdanov- 
Takens normal form 



dxi = l[x2]ds + ^^FidW^^\ 

dx2 = 7 [yi + y2X2 + xl + sxii2]ds + -^Fz^PF^^) ^ 



(80) 



with suitable slow variables y = (2/1,2/2) so that we approach the critical BT-transition at {x,y) 
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(0, 0). Consider a linear map 




where B G M^^^ is invertible. We know that 

Var(5i(y)) = O; (^^) , Var(52(y)) = (^-^) , Cov(ii(y), X2(y)) = 
as yi — 7- 0. After applying the transformation B a formal calculation yields 

Var(xi) = Var(6iixi + 612X2) = 6?iVar(xi) + 6i2Vai'(^i) + 26ii6i2Cov(xi, X2) 

Var(x2) = Var(52iii + ^22^2) = bliYa,i{xi) + 622Var(xi) + 262i622Cov(xi, X2) 

This an explanation why both variances Vi increase like 0*{l/{yi^c — yi)) in Figure [TTl The scaling 
law from the Hopf bifurcation dominates the scaling law from the saddle-node bifurcation near a 
CO dimension- two Bogdanov-Takens point when the system is not in normal form. 

We conclude this section with some potential implications for ecological modeling and ecosystem 
management. Once we have passed the BT-point the system transitions with high probability to 
a far-away equilibrium (see Figure [TOj). In particular, the density of the prey population increases 
dramatically. In this scenario it will be very difficult to reverse the system to the original state as 
the region Q2 of slow variable/parameters is very narrow near the BT-point. The most interesting 
aspect of the BT-transition in the Bazykin model (|77p is that just measuring the variances, without 
a preliminary normal form transformation, can be misleading. Measurement and fitting indicate 
a variance increase governed by C*(l/y) which could just indicate a supercritical Hopf transition 
from region Q2 to i.e. passing the (red) Hopf curve in Figure [9j This transition would not be 
critical and can easily be reversed. The slower variance increase of the critical fold transition is 
hidden near the BT-point! 
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7.5 Biomechanics and Control near Instability 




Figure 12: Panels A and B show a sketch of the Euler buckling experiment as considered in |97]. 
The force J-g compresses the spring which should stay in the upright/vertical position as shown in 
B. The bifurcation diagram on the right shows the subcritical pitchfork (|8ip with parameter values 
(j82p . The pitchfork (branch point [BP]) from the attracting equilibrium branch (think red line) 
occurs at F3 = 3.3. The unstable branches (dashed blue) undergo a further fold bifurcation (limit 
point [LP]). In A we see what happens when the spring buckles and leaves the vertical position. 

In [97] the authors investigate how humans control a spring near instability. The experimental 
setup asks participants to use their thumbs to compress the spring near the threshold of the classical 
Euler buckling instability; see Figure [T^ A mathematical model for this problem is provided 
by a subcritical pitchfork bifurcation with quintic non-linearity given by 

B' =Pi{Fs-p2)0+p^e^ -VaO"" (81) 

where J-s, Pj for j G {1, 2, 3, 4} are parameters and 9 represents the angle of the spring with respect 
to its vertical/upright position [ST]. The parameter Tg is viewed as the force applied to the spring. 
The bifurcation diagram of (j8ip is shown in Figure [T2j To stay within the framework of |97| we 
have chosen fixed parameter values 

pi = 2.639, P2 = 3.3, P3 = 106.512, p4 = 385. (82) 

The experiment in |97] asked participants to slowly compress the spring so that it does not buckle 
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Figure 13: (b)-(d) Sample paths for ([83]) with F{y) = 1 (red), F{y) = — y (green), F{y) = yc—y 
(blue) and g{x,y) = 1; fixed parameter values are given in ([52]) and {e,a) = (0.005,0.01). The 
initial condition is (a;o,yo) = (0,2). The realization of the noise W = Wg is the same for all 
three paths. In (a) we calculate the variance V = Var(x(y)) for each path using a sliding window 
technique. Note that we can already spot in the time series that variance increases for F{y) = 1, 
stays roughly constant for F{y) = yjyc — V and decays to zero for F{y) = yc — y as y tends towards 
the pitchfork critical transition at yc = 3.3. 



but also comes as close as possible to the pitchfork bifurcation. In Figure [12] this corresponds to 
moving along the stable equilibrium branch {(^,^5) G : < 3.3}. The experimental data do 
contain quite a bit of noise so that it is very reasonable to consider the system 



dx 



dy 



7 [pi{y - P2)x + P3X^ -P4X^] ds + ■^F{y)dW, 



1 ds. 



(83) 



The deterministic critical manifold of (I83p is Cq = {{x,y) £ : pi{y—p2)x-\-pzx^ — p4_x^ = 0}. We 
focus on the trivial branch Cq = {x = 0} and the attracting subset Cq := Cg H {y < 3.3}. In the 
previous applications we usually assumed that F{y) = const, which corresponds to additive noise. 
For the spring compression experiment this does not seem reasonable since participants could try 
to minimize the noisy fluctuations once they are very close the subcritical pitchfork bifurcation; in 
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Figure 14: Average variance V = Yai{x{y)) for ([83l) for F[y) = 1 (red), F{y) = ^Jyc — y (green), 
F{y) = yc — y (blue) and g{x, y) = 1 over 100 sample paths; fixed parameter values are given in 
(|82]) and (e,cr) = (0.005,0.007). The initial condition is (xo,yo) = (0,2). The point where the three 
variances cross corresponds approximately to y* = 2.3. This is expected since the three functions 
in (|84l) are equal at y = y* . 



fact, they know that a noise-induced critical transition could occur before the bifurcation point. 
Figure [T3r b)-(d) shows sample paths for different types of noise 



F{y) = 1, F{y) = Vyc - y, F{y) =y^-y. 



(84) 



We used the same realization for dW for all three paths. It can already be observed that we have 
three different behaviors ("increase, constant, decay") for the variance V = Var(x(y)). Figure [HI 
confirms this behavior as it shows the average variance over 100 sample paths for the different types 
of noise given in (|84p . We can calculate from Theorem 15.11 that to leading order in the approach 
towards the pitchfork, but not in a small neighbourhood near it, we have the scaling laws 



Var(x(y)) 



0. (EM\=o* 

^\yc-yj ^\yc-yj 



O* . 

y \y-yc 

o;{i) 



yc 



if F{y) = 1, 

if F{y) = Vyc - y, 

if F{y) =yc-y. 



This explains precisely what is shown in Figure O and shows that multiplicative noise can yield 
a wide variety of different early-warning signals or even no visible trend of the variance near a 
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critical transition. Hence we can conjecture that balancing/controlling objects near an instability 
involves suitable noisy perturbations and the quick processing of a time series history to generate 
the appropriate control. 

8 Discussion and Outlook 

This paper has only started to develop a mathematical framework for critical transitions and 
prediction. Here we briefly outline the main steps and how this framework can be extended to 
address future problems. 

The first part of this paper, motivated by Definition 12.21 only covers the singular limit e = 0, 
£7 = 0. We derive slow flow conditions to reach a critical transition and record the relevant 
linearizations to develop stochastic scaling laws. Although this is the most precise starting point 
one could consider extensions. In fact, the sample paths viewpoint of Definition 12.21 naturally 
extends. Let 

7.,<x = 7. At) ■■ [0,r] ^ M™+", 7.,<x(0) = 7(0) = (x(0),y(0)) 

be a sample path of (I27p . The first deterministic extension is to consider 70,0 but remove the 
requirement from Definition 12.21 that the transition point p is normally hyperbolic and to change 
(CI) so that a candidate 70,0(^^-1) ^j) can lie in any part of the critical manifold. This allows for 
canard orbits and delay as shown in Figure [T5lfb). As an example consider the pitchfork bifurcation 
([7]) with y(0) < then the point {x,y) = (0, min(— ?/(0), y;,)) becomes a critical transition where 
y;, > is the buffer point [76\ I77j . For delays and canards the problem of critical transitions 
becomes global in at least two ways: 

(Gl) The initial condition matters to determine which points are critical transitions. 
(G2) The global distance between critical manifolds becomes relevant. 

To understand (G2) consider the supercritical pitchfork bifurcation d?]). Depending on the 
initial condition there may be a jump at p = {xp,yp) for < yp ^ 1 or < yp = 1. The length 
of the fast segment to the next attracting critical manifold y = from p is see Figure [T5j 
Applications clearly require a case distinction between <^ 1 which is usually not viewed as 
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critical and ^Jy^ = 1 which should probably be called critical. Hence an extension to canards 
must append a global distance measure to the sample path space, e.g. the minimum or maximum 
distance from the transition point to the fast subsystem attractor; see Figure [TSTb) where canards 
with or without head usually yield two different distances. Since this paper entirely restricts to a 
local theory we do not discuss this aspect further. 




Figure 15: Illustration of possible extensions to Definition 12.21 (a) Phase space for (j27p with 
f{x,y) = yx — x^, F{x,y) = 1 and g{x,y) = 1. The critical manifold Co (grey) and two sample 
paths 7e^o and j^^a- (black) for e = 0.01 = a with initial condition (x(0),y(0)) = (0.9,-0.9) are 
shown, (b) Phase space for (p7|) with f{x,y) = y — ^ — x, g{x,y) = 1 — x and F{x,y) = 1 
with a non-generic fold at {x,y) = (1,-2/3). Again we show two sample paths (black) and the 
critical manifold (grey). Note that the path 70,0- for a = 0.5 cannot drift in y but will switch, on 
exponentially long time scales, between the two attracting branches of the critical manifold. The 
inset shows a time series for this path on a subexponential time scale. 

Another possible extension is to consider sample paths 7^,0 for < e <C 1; see Figure [TSTa). 
In this case, the extension can just be defined by requiring that (i_H'(7e,05 7o,o) — as e — 
i.e. by checking whether candidates that have a critical transition in the singular limit perturb. 
The perturbation results are known for the fold, pitchfork, transcritical and Hopf bifurcations 
[64:\ [62l [63| I76j . Partial results are available for the Bogdanov-Takens bifurcation [25j and the 
cusp [22] is work in progress; the remaining codimension-two problems are expected to be solvable 
with similar ideas. One could also add generic cases for higher-dimensional and non-minimal slow 
variables such as folded sing ularities in [90] . 

The case "jefl is primarily of mathematical interest since for 'j^^u and o" > a delay/canard effect 
is shortened substanially by noise in applications (see e.g. Theorem 2.11 of [12j) as long as the noise 
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is not exponentially small [87] . A typical delay time is of order y^e| Ino"! so the local singular limit 



results are relevant; see also Figure flST a). However, there is a major open issue for applications we 
do not address here corresponding to transitions driven purely by noise or a combination of noise 
and bifurcations. 

The case 70,0- for a sample path starting near an attracting critical manifold Cq^ is covered by 
the theory of large deviations [36j and purely noise-induced transitions can occur to an attracting 
critical manifold Cq^; see Figure [TSlb) where the upper path will eventually escape. Again, the 
sample path viewpoint is well-suited as we ask for an estimate of probabilites e.g. 



for suitable small constants 81^2 and a given time t* > 0. Hence one can again use paths and 
Definition 12.21 as a basis but then has to add for each point on Cq ^ a probabilistic description how 
likely the escape is which usually yields exponentially long time scales to escape. This is again a 
global problem. For cases with one fast variable and e = it is often possible to obtain explicit 
solutions using Fokker-Planck equations e.g. see [37 } 111 IT U l9T | l65]. 

Remark: After the suggestion of the Definition 12.21 in [65], recent work of Ashwin et al. [6] 
suggested a related applied classification of critical transitions distinguishing between B-tipping 
('bifurcation-induced'), N-tipping ('noise- induced') and R-tipping ('rate- induced'). Basically B- 
tipping aims to cover paths 7^,0 for e — )• and N-tipping considers paths 70,0-; it is currently work 
in progress to understand R-tipping better. 

The most general case is to consider j^^a- for o", e > where noise-induced escape shortly before 
a fast subsystem bifurcation point on non- exponential time scales becomes relevant. One of the 
key goals of the mathematical framework presented in this paper was to also allow for a natural 
extension of the methods and definitions to this case. It is future work to combine the ideas from 
Definition 12 . 21 by adding to it pathwise probability estimates of the form (IBSp . This should yield the 
full mathematical framework based upon sample paths with all parameters: o" > 0, e > 0, distance 
to the next attractor and escape probability during [0,T]. 



P inf t : dh{joAt),CSn < 52,rf//(7o,a(0),Co"^) < <5i >f 




) 



(85) 
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For the local theory of codimension-one bifurcations several studies on various regimes with 
(7, e > near bifurcation points exist. Overall, the fold [iniET], pitchfork/transcritical [Ml US] and 
Hopf bifurcations [HI [891 [16] quite well understood. One basic insight is that scaling regimes are 
identified under which noise-induced effects or deteterministic drift dominate. Another important 
conclusion are probabilistic estimates for certain distinct dynamical regimes to occur. For higher 
codimension phenomena not many results are known but see e.g. \n\ I18|. As far as the stochastic 
scaling laws for codimension- two cases considered in this paper are concerned there does not seem 
to be any work prior to this paper in this direction. 

The second contribution of this paper is to understand fluctuations and scaling laws of paths 
better before fast subsystem bifurcations to determine early-warning signs. In particular, leading- 
order scaling behaviour for covariance matrices have been derived. We have only covered the basic 
case of local bifurcations up to codimension- two with white noise in the region (Rl) with a suitable 
scaling of noise and time scale separation which makes early escapes unlikely. Large fluctuations 
before the bifurcation and scaling results near bifurcations are certainly not well-studied for all 
bifurcations up to codimension two. Early-warning signs for other types of noise (colored noise, 
shot/burst noise [37]), for degenerate noise terms [92j and for more general stochastic processes 
(e.g. Levy Processes [581 [52]) interesting directions. As before, sample paths and singular limits 
are still available, even for very general high-dimension bifurcations and stochastic processes. 

Global bifurcations [69] have not been considered and would be an interesting direction for 
future analysis. There is work in progress to understand these bifurcations and their warning signs 
in models as well as in a normal form setup. Another possible extension are early-warning signs 
for spatially-extended problems; see \29\ [321 |27] for models from ecology. In this context, it is 
well-known that many classes of pattern-forming partial differential equations (PDEs) and stochas- 
tic partial differential equations (SPDEs) can be written as evolution equations with well-defined 
paths or stochastic sample paths [48[ I81|. Several relevant PDEs, such as excitable systems [75| \7\ 
with diffusion, are often already in a natural fast-slow form. Presumably one should find many 
other interesting early-warning signs for spatial systems but these could also be more difficult to 
measure and apply in practical applications since the collection and analysis of much larger data 
sets arises; a typical area where this already proved to be very difficult are epileptic seizures [7^1172]. 
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The third contribution of the current work are examples, several of them in application domains 
(epidemics, systems biology and biomechanics) where the new techniques for early-warning signs 
have not been considered. Furthermore, the examples provide illustrations of the theory and also 
show its limitations where prediction becomes impossible or misleading if one relies on the scaling of 
the variance. There are many important directions for making the theory more applicable e.g. de- 
tailed statistical tests such as receiver-operator curves [Ml [671 [IS] i analysis of limited data and its 
interpretation [31], linking critical transitions to experiments [6'6\ I98j. desirable tipping points in 
applications and their control [551 [66] as well as networks and deterministic metastability |67j . 
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